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ABSTRACT 


A  Navier-Stokes  problem  solver,  developed  by  L.  N.  Sankar,  is  installed  and 
verified  on  the  NASA  Ames  Cray  X/MP-48  computer  and  is  used  to  calculate  the  flow 
field  about  a  NACA  0012  airfoil  oscillating  in  pitch.  Surface  pressure  distributions  and 
integrated  lift,  pitching  moment,  and  drag  coefficients  versus  angle  of  attack  are 
compared  to  existing  experimental  data  for  two  cases,  involving  deep  dynamic  stall  and 
fully  attached  flow  at  and  below  a  freestream  Mach  number  of  .3.  The  flow  field  about 
the  oscillating  airfoil  is  investigated  through  the  study  of  contour  plots  of  pressure, 
density,  Mach  number,  and  stream  function.  The  effect  of  turbulence  modeling  is 
explored  through  use  of  the  Baldwin-Lomax  model  and  a  modification  designed  to 
prevent  underprediction  of  maximum  lift.  Finally,  Reynolds  number  and 
compressibility  effects  are  investigated  by  repeating  the  deep  stall  simulation  at  one- 
tenth  the  experimental  Reynolds  number  and  Mach  numbers  of  .3  and  .5.  The  latter 
conditions  are  intended  for  comparison  with  the  results  of  wind  tunnel  experiments 
being  planned  at  NASA  Ames  Fluid  Mechanics  Laboratory. 


4 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . . 

II.  MATHEMATICAL  FORMULATION . 17 

A.  THE  SOLUTION  DOMAIN . . 

1.  Algebraic  Grid  Generation  . 17 

2.  Initial  Conditions  and  Boundary  Conditions . 18 

B.  GOVERNING  EQUATIONS . 23 

C.  NUMERICAL  SOLUTION  . 26 

1.  Linearization  and  Factorization . 27 

2.  Application  of  Boundary  Conditions  . 29 

3.  Artificial  Dissipation . . 

D.  TURBULENCE  MODEL  . . 

1.  Baldwin- Lomax  Model  . . 

2.  Modified  Model . . . 33 

III.  DESCRIPTION  OF  THE  CODE . 34 

A.  MAIN  PROGRAM  . . 

B.  GRID  GENERATION  . . 

C.  COMPUTATIONAL  STEPS . . 

IV.  RESULTS  AND  DISCUSSION . . 

A.  STEADY-STATE  CALCULATIONS . 44 

B-  COMPARISON  WITH  DEEP  STALL  EXPERIMENTAL 

L/A  1 A  . . . . 

c  flow . 67 

CONATIONS5  UNDER  PROSED  EXPERIMENTAL  ^ 

V.  CONCLUDING  REMARKS . 80 

APPENDIX  A:  SANKAR  NAVIER-STOKES  SOLVER . 82 


5 


APPENDIX  B:  NOTES  ON  USE  OF  THE  NAVIER-STOKES  SOLVER . 113 

1.  JOB  CARDS . 113 

2.  MAIN  PROGRAM  . 114 

3.  DATACARDS . 114 

4.  ADDITIONAL  NOTES . 116 

LIST  OF  REFERENCES . 117 

INITIAL  DISTRIBUTION  LIST . 119 


6 


LIST  OF  TABLES 


1.  SUBROUTINE  CALLS  FROM  THE  MAIN  PROGRAM  . 35 

2.  GRID  GENERATION  SUBROUTINES  CALLED  . 37 

3.  SUBROUTINES  USED  IN  COMPUTATIONAL  STEPS  . 39 

4.  INPUT  CONDITIONS  FOR  COMPUTER  SIMULATIONS  . 43 


7 


LIST  OF  FIGURES 


1.1  Dynamic  Stall  Chronology  (from  Carr,  Ref.  1) . 15 

2.1  Symmetric  Airfoil  in  the  Physical  Plane  (top)  Unwrapped  in  the 

Intermediate  Plane  (bottom) . 19 

2.2  Construction  of  the  Grid  in  the  Intermediate  Plane  (top) 

Corresponding  Computational  Plane  Grid  (bottom) . 20 

2.3  Algebraic  C-Grid  in  the  Physical  Plane  . 21 

2.4  Detail  of  Grid  near  the  Airfoil  . 21 

4.1  Steady-State  Attached  Flow  Mao*  .3,  Re*  LOO  x  106,  a*  12°  Top- 

Density  (1),  Pressure  (r)  Bottom-Mach  No.  (1),  Stream  Function  (r) . 45 

4.2  Steady-State  Separated  Flow  Mqo *  .3,  Re*  1.00 x  106,  a*K° 

Density  (1)  and  Pressure  (r)  at  5100,  5610,  and  6120  time  steps . 46 

4.3  Steady-State  Pressure  Distribution  at  Maximum  Lift  Mqq  *  .301, 

Re*  3.91  x  106,  a  -  13.5°  Inset-Experimental  Data  (Ref.  8) . 47 

4.4  Lift  and  Pitching  Moment  Coefficients-Baldwin- Lomax  Model 

Moo *.283,  Re* 3.45  x  106,  a*  15°— 10°cos(.3t) . 49 

4.5  Pressure  Drag  Coefficient-Baldwin-Lomax  Model  Mqq  *  .283, 

Re*  3.45  x  106,  a*  15B-10Bcos(.3t) . 50 

4.6  Surface  Pressure  Coefficient  Carpet  Plots-Baldwin-Lomax  Model 

Moo*. 283,  Re* 3.45  x  I06,  a*  15#-10°cos(.3t),  Top-theoretical, 
Bottom-experimental  (Ref.  8)  . 51 

4.7  Density  Contour  Plots,  Upstroke-Baldwin-Lomax  Model  M oo  *  .283, 

Re-  3.45  x  106,  a-  15°-10°cos(.3t)  6.31,  9.95,  14.93,  19.99,  23.65,  25 
degrees  .  52 

4.8  Density  Contour  Plots,  Downstroke-Baldwin-Lomax  Model 

Mqo -.283,  Re- 3.45  x  106,  a  -  15°-10°cos(.3t)  23.67,  20.02,  15.03, 

10.03,  6.36,  5  degrees . 53 

4.9  Pressure  Contour  Plots,  Upstroke-Baldwin-Lomax  Model 

Moo -.283,  Re- 3.45  x  106,  a*  15°-10°cos(.3t)  6.31,  9.95,  14.93, 

19.99,  23.65,  25  degrees . 54 


8 


4.10  Pressure  Contour  Plots,  Downstroke- Baldwin- Lomax  Model 

Mqo -.283,  Re- 3.45  *  106,  a  -  15°-10°cos(.3t)  23.67,  20.02,  15.03, 

10.03,  6.36,  5  degrees . .  55 

4.1 1  Mach  Number  Contour  Plots,  Upstroke-Baldwin- Lomax  Model 
Moo -.283,  Re- 3.45  x  106,  a- I5°-10°cos(.3t)  6.31,  9.95,  14.93, 

19.99,  23.65,  25  degrees . 55 

4.12  Mach  Number  Contour  Plots,  Downstroke--BaIdwin-Lomax  Model 
Moo  -  .283,  Re  -  3.45  x  106,  a  -  15°-10°cos(.3t)  23.67,  20.02,  15.03, 

10.03,  6.36,  5  degrees . . 

4.13  Stream  Function  Plots,  Upstroke-Baldwin-Lomax  Model  Mqo  -  .283, 

Re-  3.45  x  io6,  a  -  15°-10°cos(.3t)  6.31,  9  95,  14.93,  19.99,  23.65,  25 
degrees .  5g 

4.14  Stream  Function  Plots,  Downstroke-Baldwin-Lomax  Model 

Moo -.283,  Re- 3.45  x  106,  a  -  15°-10°cos(.3t)  23.67,  20.02,  15.03, 

10.03,  6.36,  5  degrees . . 

4.15  Lift  and  Pitching  Moment  CoelFicients-Modified  Model  until  25° 

Mqo -.283,  Re- 3.45  x  106,  a- 15°-10°cos(.3t) . 62 

4.16  Surface  Pressure  Coefficient  Carpet  Plots-Modified  Model  until  25° 

M» *-283,  Re- 3.45  x  106,  a  -  15°-- 10°cos(.3t),  Top-theoretical, 
Bottom-experimental . . 

4.17  Density  (top)  and  Pressure  (bottom)  Plots-Modified  Model  until  25° 

Moo -  283,  Re- 3.45  x  IO6,  a-  15°-10°cos(.3t)  25°  (1),  23.67°  (r) . 64 

4.18  Mach  Number  and  Stream  Function  Plots-Modified  Model  until  25° 

Moo  -  -283.  Re-  3.45  x  106,  a  -  15°-10°cos(.3t)  25°  (1),  23.67°  (r) . 65 

4.19  Lift  and  Pitching  Moment  Coefficients-Modified  Model  until  19° 

Mqo  -  -283,  Re-  3.45  x  106,  a  -  15°— 10°cos(.3t) . 66 

4.20  Lift  and  Pitching  Moment  Coefficients-Baldwin-Lomax  Model 

Moo-184,  Re-2.45  x  106,  a- 8°— 10°cos(.4t) . 68 

4.21  Surface  Pressure  Coefficient  Carpet  Plot- Baldwin- Lomax  Model 

Moo 35 .184,  Re- 2.45  x  106,  ct  —  8°— 10°cos(.4t)  Top— theoretical, 
Bottom-experimental  (Ref.  8)  . . 

4.22  Flow-field  Plots-Baldwin-Lomax  Model  Mqo  =  .184,  Re  -  2.45  x  106, 
a  -  8°-10°cos(.4t)  Top-Density  (1),  Pressure  (r)  Bottom-Mach  No. 

(1),  Stream  Function  (r) . . 


9 


4.23  Lift  and  Pitching  Moment  Coefficients-Modified  Turbulence  Model 


Mqo  =  .184,  Re=2.45x  106,  a»  8°-10°cos(.4t) . 72 

4.24  Surface  Pressure  Coefficient  Carpet  Plots- Modified  Turbulence 
Model  Mx  =  .184,  Re=  2.45  x  106,  a  =  8°-10°cos(.4t)  Top- 

theoretical, Bottom-experimental  (Ref.  8)  . 73 

4.25  Flow-field  Plots-Modified  Turbulence  Model  Mqo  =  •  184, 

Re=  2.45  x  106,  a  =  8°-10#cos(.4t)  Top-Density  (1),  Pressure  (r) 

Bottom-Mach  No.  (1),  Stream  Function  (r)  . 74 

4.26  Lift  and  Pitching  Moment  Coeflficients-Baldwin-Lomax  Model 

Moo  =  284,  Re=.345x  106,  a=  15°-10°cos(.3t) . 76 

4.27  Surface  Pressure  Coefficient  Carpet  Plots-B^ldwin-Lomax  Model 
Top:  Moo  =  -284,  Re  =  .345  x  106,  a»  15°-I0°cos(.3t)  Bottom: 

Re=  3.45  x  106(Previous  Results) . 77 

4.28  Lift  and  Pitching  Moment  Coefiicients-Baldwin-Lomax  Model 

Moo  =  .5,  Re  =  .345  x  106,  <*=  15°-10°cos(.3t) . 78 

4.29  Surface  Pressure  Coefficient  Carpet  Plots-Baldwin-Lomax  Model 
Top:  Moo  =  .5,  Re=.345  x  106,  a=  15#-10°cos(.3t)  Bottom: 

M  oo  =  -284  (Previous  Results)  . 78 


10 


TABLE  OF  SYMBOLS 


A,  B 

a 

c 

Cd 

q 

Cm 


D 

e 

F,G 

F*.  G* 

I 

J 

k 

C 

Moo 

Pr 


P 

q 

q* 

QX.  Qy 
R,  s 
R*,  S* 
Re 
U,  V 
u,  v 
X,  y,  t 
z 


a 

a 

a 


o 

1 


Jacobian  matrices 
speed  of  sound 
chord  length 
drag  coefficient 
lift  coefficient 

pitching  moment  coefficient 

pressure  coefficient 

Van  Driest  damping  function 

total  energy  per  unit  volume 

flux  vectors  in  the  physical  plane 

flux  vectors  in  the  transformed  plane 

identity  matrix 

transformation  Jacobian 

reduced  frequency 

eddy  viscosity  mixing  length 

freestream  Mach  number 

Prandtl  number 

pressure 

unknown  flow-field  vector  in  the  physical  plane 

unknown  flow-field  vector  in  the  transformed  plane 

heat-flux  vector  c-.  aponents 

viscous  stress  vectors  in  the  physical  plane 

viscous  stre: .  vectors  in  the  transformed  plane 

Reynolds  number 

contravariant  velocity  components 

velocity  components  in  the  physical  plane 

coordinates  in  the  physical  plane 

complex  coordinate  representation  of  a  point  in  the  physical  plane 

angle  of  incidence 

the  mean  angle  of  incidence 

the  amplitude  of  oscillations 


11 


Y 

specific  heat  ratio 

eE’eI 

coefficients  of  explicit  and  implicit  artificial  dissipation  terms 

P 

viscosity 

eddy  viscosity 

5.  n.  t 

coordinates  in  the  transformed  plane 

nx.  • 

transformation  metrics 

p 

density 

*xx’  *xy’  *yy 

components  of  shear-stress  tensor 

(0 

vorticity  or  circular  frequency 

Operators: 

d 

6 

\ 

V 


partial  derivative 

central  difference 

incremental  change 

backward  difference  or  del  operator 


12 


ACKNOWLEDGEMENTS 


The  research  for  this  thesis  was  conducted  at  the  Fluid  Mechanics  Laboratory  of 
the  NASA  Ames  Research  Center  under  the  Navy-NASA  Joint  Institute  of 
Aeronautics.  This  work  is  related  to  the  research  program  being  carried  out  for  the 
Army  Research  Office  under  proposal  number  23394-EG.  Work  was  performed  under 
the  guidance  of  Professor  Satyanarayana  Bodapati,  thesis  advisor,  and  Dr.  Lawrence 
W.  Carr  of  NASA  Ames,  co-advisor.  I  would  like  to  express  my  sincere  gratitude  to 
them  for  providing  me  with  the  opportunity  to  become  acquainted  with  the  field  of 
computational  fluid  dynamics  and  for  their  patience  and  encouragement. 

I  would  also  like  to  thank  Professor  L.  N.  Sankar  of  the  Georgia  Institute  of 
Technology  for  providing  the  computer  program  used  and  for  being  so  generous  with 
his  time  during  several  telephone  conversations  and,  along  with  his  associate  Dr.  Wei 
Tang,  during  a  personal  visit. 

My  sincerest  appreciation  goes  to  Joan  Thompson  of  Sterling  Software  for  the 
professional  and  friendly  manner  in  which  she  assisted  me  on  occasions  too  numerous 
to  recall  in  working  with  the  Ames  computer  facilities;  and  to  Rosalie  Lefltowitz,  also 
of  Sterling  Software,  for  her  help  with  graphical  presentation,  which  included  providing 
the  carpet  plotting  routine  used  for  this  study. 

Finally,  I  wish  to  thank  my  wife,  Barbara,  for  both  her  patience  and  her 
invaluable  help  in  preparation  of  this  thesis,  and  my  daughter,  Andrea,  for  giving  up 
her  father  during  the  final  few  months  of  this  study. 


13 


I.  INTRODUCTION 


Dynamic  stall  refers  to  the  delay  of  stall  onset  beyond  static  stall  angles 
experienced  by  airfoils  and  wings  in  unsteady  motion.  The  phenomenon  first  drew 
serious  attention  in  connection  with  helicopter  aerodynamics  when  conventional 
methods  of  analysis  proved  incapable  of  accurately  predicting  performance  for  vehicles 
in  high  speed  forward  flight.  The  observed  increase  in  overall  lift  could  be  explained  if 
the  lift  on  the  blade  moving  opposite  to  flight  direction  was  greater  than  predicted  by 
steady  flow  calculations.  [Ref.  1)  It  was  experimentally  observed  that  the  lifting 
characteristics  of  rotors  could  be  adequately  modeled  by  an  airfoil  oscillating  in  pitch. 
The  basic  mechanism  is  the  shedding  of  a  strong  leading-edge  vortex  which  distorts  the 
pressure  distribution  as  it  travels  over  the  upper  surface  of  the  airfoil  and  leads  to  the 
abrupt  changes  in  lift  and  pitching  moment  usually  associated  with  dynamic  stall. 
[Ref.  2) 

In  addition  to  helicopter  rotors,  the  dynamic  stall  phenomenon  is  observed  in 
such  diverse  aerodynamic  applications  as  wind  turbines,  jet  engines,  and  rapidly 
maneuvering  aircraft.  Interest  in  expanding  the  design  envelopes  of  fighter  aircraft 
through  both  post  stall  maneuvers  and  extended  conventional  maneuverability  has 
increased  interest  in  dynamic  stall  research.  Progress  in  this  area  depends  upon 
improved  knowledge  of  details  of  viscous  flow  over  airfoils  in  dynamic  stall.  [Ref.  1] 
The  basic  airfoil  dynamic  stall  process  is  depicted  in  Figure  1.1,  taken  from  Reference 

1.  The  flow-field  pattern  at  any  instant  is  critically  dependent  upon  the  entire  time 
history,  so  understanding  the  chronology  of  events  including  time  prior  to  the  vortex 
shedding  is  crucial.  There  are  four  basic  stages; 

1.  trailing-edge  flow  reversal,  progressing  forward  along  the  airfoil, 

2.  formation,  growth,  and  shedding  of  the  vortex, 

3.  full  stall, 

4.  boundary  layer  reattachment. 

The  specific  characteristics  of  a  given  oscillating  airfoil  depend  primarily  on 

1.  Mach  number, 

2.  Reynolds  number, 

3.  reduced  frequencv,  k.  the  ratio  of  the  vertical  velocity  of  the  leading  edge  to  the 
freestream  velocitv,  k  =  coc/2U,  where  to  =  circular  frequency,  c  =  airfoil 
chord,  and  U  =  freestream  velocity, 
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PITCHING  MOMENT.  CM  NORMAL  FORCE. 


0  5  10  16  20  25 

INCIDENCE,  a,  d*g 


(•)  STATIC  STALL  ANGLE  EXCEEDED 


(b)  FIRST  APPEARANCE  OP  FLOW 
REVERSAL  ON  SURFACE 


(c)  LARGE  EDDIES  APPEAR  IN 
BOUNDARY  LAYER 


(d)  FLOW  REVERSAL  SPREADS  OVER 
MUCH  OF  AIRFOIL  CHORD 


(•)  VORTEX  FORMS  NEAR 
LEADING  EDGE 


(f)  LIFT  SLOPE  INCREASES 


(9)  MOMENT  STALL  OCCURS 


(W  LIFT  STALL  BEGINS 


(I)  MAXIMUM  NEGATIVE  MOMENT 


(j)  FULL  STALL 


00  BOUNDARY  LAYER  REATTACHES 
FRONT  TO  REAR 


(I)  RETURN  TO  UNSTALLED  VALUES 


Figure  1.1  Dynamic  Stall  Chronology  (from  Carr,  Ref.  1). 
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4.  mean  angle  of  attack,  aQ, 

5.  oscillation  amplitude,  ctj.  [Refs.  1,3] 

While  much  of  previous  research  has  involved  experiment;,,  primarily  on  airfoils 
oscillating  in  pitch,  advances  in  computer  technology  have  increased  the  potential  of 
computational  methods  to  model  details  of  unsteady  flow  fields  accurately.  Early 
analytic  approaches  based  on  inviscid  models  or  semi-empirical  analysis  were  very 
limited  due  to  the  complexity  of  the  dynamic  stall  phenomenon  and  its  interrelation 
with  viscous  effects.  Moreover,  boundary  layer  corrections  are  not  appropriate  due  to 
the  presence  of  large  regions  of  separation.  [Ref.  4]  The  Navier-Stokes  equations  can 
deal  with  flow  separation  and  shock-boundary  layer  interaction,  and  the  thin-layer 
approximation  has  been  successfully  applied  to  airfoils  near  maximum  lift,  including- 
though  somew'hat  less  successfully-conditions  beyond  maximum  lift  coefficient  [Ref.  5]. 
The  thin-layer  Navier-Stokes  equations  may  not  be  appropriate  in  the  case  of  highly 
separated  flow's,  howrever,  and  in  the  case  of  dynamic  stall,  the  viscous  layer  is  of  the 
same  order  as  the  airfoil  chord  during  the  shedding  of  the  strong  leading-edge  vortex. 
The  full  Navier-Stokes  equations  must  therefore  be  solved.  [Ref.  2] 

A  Navier-Stokes  problem  solver  has  been  developed  by  L.  N.  Sankar  and  his 
associates  a*,  the  Georgia  Institute  of  Technology  and  was  made  available  for  the 
present  study.  An  implicit  finite-difference  procedure  is  used  to  solve  the  two- 
dimensional,  Reynolds-averaged,  compressible  Navier-Stokes  equations  in  strong 
conservation  form.  The  flow  solver  features  a  moving,  body-fitted  coordinate  system, 
an  algebraic  turbulence  model  (Baldwin-Lomax),  and  an  alternating-direction-implicit 
(ADI)  time-marching  algorithm.  It  can  also  be  u?ed  in  an  inviscid  mode  to  solve  the 
Euler  equations.  [Refs.  6,7]  The  goals  of  the  present  study  were: 

1.  Develop  procedures  for  exercising  the  code  on  the  Cray  X/MP-48  computer  .to 
obtain  output  in  a  form  suitable  for  comparing  the  pressure  distribution  with 
experimental  results  and  for  investigating  details  of  the  computed  flow  field. 

2.  Apply  the  code  to  steady-state  cases  and  to  the  conditions  of  existing  dynamic 
stall  experimental  data  [Ref.  8]  and  investigate  the  effect  of  varying  input 
parameters. 

3.  Obtain  flow-field  solutions  under  proposed  conditions  of  a  series  of  wind  tunnel 
experiments  to  be  conducted  at  the  Fluid  Mechanics  Laboratory  of  the 
National  Aeronautics  and  Space  Administration,  Ames  Research  Center. 
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II.  MATHEMATICAL  FORMULATION 
A.  THE  SOLUTION  DOMAIN 

The  first  step  in  any  numerical  solution  procedure  is  to  discretize  the  physical 
domain.  This  involves  the  selection  of  suitable  boundaries  and  the  choice  of 
appropriate  nodal  points  within  the  boundaries  to  describe  the  characteristics  of  the 
physical  domain  accurately.  In  the  solver  used  here  a  solution  grid  is  constructed 
about  the  airfoil,  which  is  then  rotated  along  with  the  airfoil  for  dynamic  calculations. 
With  the  solution  domain  defined,  it  is  then  necessary  to  specify  conditions  at  the 
boundaries.  For  time  dependent  problems  such  as  those  considered  here,  initial 
conditions  at  the  nodal  points  must  also  be  defined. 

1 .  Algebraic  Grid  Generation 

A  mathematical  transformation  from  the  physical  plane  to  a  computational 
plane  introduces  significant  advantages  for  numerical  solution.  Chief  among  these  is 
the  fact  that  the  physical  boundaries  are  mapped  into  rectangular  surfaces  in  the 
transformed  plane.  Grid  points  can  be  concentrated  in  physical  regions  experiencing 
the  largest  gradients  while  maintaining  the  simplicity  of  uniform  spacing  in  the 
computational  plane.  [Refs.  9,10:  pp.  519-520]  Proper  design  of  the  grid  is  crucial  to  a 
stable,  accurate  solution.  The  version  of  Sankar's  program  used  for  this  study  employs 
an  algebraically-generated  C-grid.  This  grid  is  quickly  produced,  easily  rotated  by 
simple  sine/cosine  relationships  with  the  angle  of  attack,  and  allows  a  high  degree  of 
clustering  in  the  normal  direction  to  cover  adequately  the  thin  boundary  layer  of  high 
Reynolds  number  flows. 

The  procedure  used  generates  a  sheared  parabolic  coordinate  system  based  on 
an  input  airfoil  shape.  Two  points  are  first  found  on  the  airfoil:  point  N  halfway 
between  the  nose  and  the  center  of  curvature  of  the  nose,  and  point  T  at  the  trailing 
edge.  The  cut  along  the  airfoil  wake  is  chosen  to  be  tangent  to  the  mean  camber  line 
at  the  trailing  edge.  The  airfoil  shape  and  wake  are  then  unwrapped  to  form  a  surface 
S(X)  in  an  intermediate  X-Y  plane  using  the  following  transformation: 


X  +  iS(X)  = 


(eqn  2.1) 
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where  z  =  x  +  iy  locates  a  point  in  the  physical  (x-y)  plane.  This  surface  will  then  lie 
along  the  X-axis  with  a  shallow  bump  at  the  origin  which  coincides  with  the  singular 
point  Z\^.  Figure  2.1  illustrates  this  mapping  for  a  symmetric  airfoil.  Cubic 
interpolation  is  used  to  identify  additional  points  and  smooth  the  surface  S(X).  Lines 
parallel  to  the  Y-axis  and  lines  equidistant  from  the  surface  S(X)  are  then  constructed 
before  the  sheared  cartesian  grid  is  mapped  back  to  the  physical  plane.  Each  point  in 
the  physical  domain  is  assigned  to  a  corresponding  point  in  the  computational  (vH) 
plane  through  the  intermediate  plane  relationship.  The  intermediate  and 
computational  planes  are  illustrated  in  Figure  2.2.  Unit  spacing  is  assumed  in  the 
computational  plane,  simplifying  calculations.  For  solutions  of  viscous  flows,  the  ^  = 
constant  lines  are  retained  w’hile  the  q  =  constant  lines  in  the  physical  plane  are 
redistributed  to  place  points  within  the  boundary  layer.  The  first  point  off  the  solid 
surface  is  placed  at  a  distance  specified  by  the  user  and  the  remaining  lines  are 
exponentially  distributed  to  the  far-field  boundary.  [Refs.  6,7:  pp.  40-44] 

The  final  mesh  used  is  shown  in  Figure  2.3.  The  number  of  points  defined  in 
the  vdirection  is  161  of  which  159  are  used  for  calculations,  and  41  points  are  defined 
in  the  q-direction.  A  detail  of  the  grid  around  the  airfoil  is  showm  in  Figure  2.4.  The 
spacing  along  constant-^  lines  is  relatively  fine  near  the  leading  edge,  where  the 
sharpest  gradients  are  expected,  and  coarse  at  the  trailing  edge.  In  order  to  perform 
calculations  in  the  computational  plane,  a  transformation  is  required  to  map  each 
point  to  the  corresponding  point  in  the  physical  plane.  In  the  present  solver,  a 
numerical  approach  is  taken:  finite  differences  are  used  to  relate  the  mesh  spacings  in 
the  two  planes  through  transformation  metrics.  This  method  is  outlined  in  later 
sections. 

2.  Initial  Conditions  and  Boundary  Conditions 

The  initial  conditions  for  steady-state  calculations  are  assumed  to  be  the 
freestream  conditions,  and  inaccuracies  introduced  by  this  crude  starting  approximation 
are  removed  by  viscous  effects  after  the  airfoil  is  impulsively  started.  (Inviscid 
calculations  require  proper  boundary  conditions  and  artificial  dissipation  to  correct  the 
solution).  [Ref.  7:  p.  13]  In  the  case  of  dynamic  calculations,  the  initial  conditions  are 
obtained  from  an  asymptotically  converged  steady-state  solution  at  the  minimum 
airfoil  angle  of  incidence.  The  solution  is  then  marched  in  time  as  the  oscillation  cycle 
begins  from  this  angle.  Boundary  conditions  are  explicitly  applied  at  each  time  step  on 
the  solid  boundary,  at  the  far-field  boundary,  and  in  the  airfoil  wake,  where  the  grid 
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Figure  2.1  Symmetric  Airfoil  in  the  Physical  Plane  (top) 
Unwrapped  in  the  Intermediate  Plane  (bottom). 
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generation  procedure  introduced  a  "cut"  in  the  physical  domain  between  the  airfoil 
trailing  edge  and  the  downstream  boundary.  (See  Fig.  2.3)  [Refs.  6,11] 

B.  GOVERNING  EQUATIONS 

The  unsteady,  compressible,  two-dimensional  Navier-Stok.es  equations  are  written 
in  cartesian  coordinates  as: 


dq 

IT  +  *  •  p  -  o. 


pu 

1  + 

pv 

Pu2  +  P-'„ 

PUV-'xy 

P“V-txy 

P  +  p-t„ 

(e  +  p)u.ut„.vr<y  +  Qx 

(=+P)v-vt„-ut<y+Qy 

j> 


pu 

pv 


(eqn  2.2) 


where  p,  u,  v,  and  e  are  density,  velocity  components,  and  total  energy  per  unit 
volume;  txx,  tyy,  and  txy  are  components  of  the  shear-stress  tensor  derived  by 
application  of  Stokes'  hypothesis;  and  Qx  and  Qy  are  the  heat-flux  vector  components 
given  by: 


Q.  -  — kd.T  - - — &(a2) 

1  1  (y-l)Pr  *'  ' 


(eqn  2.3) 


in  which  p  is  the  coefficient  of  viscosity,  ?r  is  the  Prandtl  number,  and  a  is  the  speed  of 
sound.  [Refs.  12,10:  pp.  480-481]  The  specific  heat  ratio,  y,  is  assumed  to  be  1.4. 
Expanding  and  rearranging  Equation  2.2  yields  [Refs.  6,7:  pp.  9-10]: 


dtq  +  axF  +  5yG  -  3XR  +  dyS 


(eqn  2.4) 
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with 
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and 


R4  "  uTxx  +  VTxy  +  K5x<a2) 
S4  *  utxy  +  Vtyy  +  Rdy(a2) 


(y-l)Pr 


(eqn  2.5) 


In  this  form  F  and  R  contain  the  flux  and  viscous  stress  terms  in  the  x-direction  and  G 
and  S  the  flux  and  viscous  stress  terms  in  the  y-direction.  Note  that  the  viscous  terms 
are  contained  on  the  right-hand  side  (setting  them  equal  to  zero  produces  the  Euler 
equations).  The  governing  equations  are  given  in  strong  conservation  form.  That  is, 
the  coefficients  of  the  derivative  terms  are  constant,  and  the  equations  may  be  written 
in  a  form  expressing  the  divergence  of  the  physical  quantities  of  mass,  momentum,  and 
energy  (Eqn.  2.2).  Such  equations  may  be  expressed  in  a  corresponding  integral  form 
through  the  divergence  theorem: 

JvV  •  Pdv  -  0  -  J$P  *  nds  (eqn  2.6) 

The  differential  statement  of  the  governing  equations  is  not  valid  across  discontinuities 
(ie.,  shocks);  however  the  integral  expression  is.  [Ref.  13:  pp.  406-408]  Using  the 
conservative  form  means  that  the  finite-difference  approximation  will  ensure 
conservation  of  the  physical  properties  and  thus  satisfy  the  corresponding  integral  form 
of  the  governing  equations.  [Refs.  14,10:  pp.  50-52]  This  gives  weak  solutions  (solutions 
of  partial  differential  equations  that  include  discontinuities)  [Ref.  10:  pp.  139-141]  in 
the  vicinity  of  shocks.  Such  methods  are  called  shock  capturing,  and  their  use  is 
important  in  treatment  of  the  compressible  dynamic  stall  problem. 
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As  discussed  earlier,  considerable  simplification  of  the  numerical  solution  may  be 
realized  by  employing  a  transformation  into  a  rectangular  plane.  It  can  be  shown  that 
the  strong  conservation-law  form  may  be  restored  following  such  a  transformation 
[Refs.  14,9:  p.  254].  The  general  transformation  is  given  by:  [Ref.  6] 


\  "  4(x,y,t) 

n  -  n(x.y.t) 

T  *  t  (eqn  2.7) 


and  the  Jacobian  and  metrics  of  the  transformation  by: 


3  *  Sxny-nx5y  -  i/(x$yn-xny£) 


(eqn  2.8) 


■  -Jxn 

*»t  ”  “XT^x""yT^y 

n  = 
ny  -  Jx^ 

*it  =  -*x\  -  yT  ny 


(eqn  2.9) 


Here  the  subscripts  indicate  partial  differentiation  with  respect  to  the  subscripted 
variable.  The  Jacobian  represents  the  magnification  factor  of  mesh  elements  between 
the  physical  and  computational  planes  [Ref.  10:  p.  530].  The  time-  and  spatial- 
derivatives  of  any  quantity  q  are  given  by: 


9t  *  +  +  Vt 

%  “  +  V* 


(eqn  2.10) 


The  governing  equation  (Eqn.  2.4)  in  the  transformed  plane  may  then  be  expressed  as: 


dxq*  +  d%¥*  +  G*  -  ^R*  +  ^S* 


(eqn  2.11) 
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where 


q*  -  q/J 

F*  -  ($XF  +  $yG  +  $tq)/J 
g*  -  (nxF  +  nyc  +  ntq)/J 
R*  -  (£XR  +  ^ys)/j 

S*  =  (nxR  +  »lyS)/J  (eqn  2.12) 

This  constitutes  a  coupled,  highly  non-linear  system  of  equations.  The  flux  vectors  F* 
and  G*  may  be  stated  in  a  less  complex  form  by  introducing  the  contravariant  velocity, 
which  is  required  to  compensate  for  grid  motion.  The  contravariant  velocity 
components  are  given  by  [Ref.  9]: 

U  *  5t  +  Sxu  +  V  -  4x(u-xt)  +  5y(v-yt) 

v  "  nt  +  nxu  +  nyv  *  nx(u-xt)  +  ny(v-yT)  (eqn  2.13) 

Using  these  velocities,  F*  and  G*  may  be  rewritten,  as  used  in  the  present 
implementation: 


1 
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l 

G*  -  — 
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-pv 

puU  +  5xp 

puv  +nxP 

pvU  +  $yp 

pvv  +qyp 

_(e  +  p)U-ntP 

_(e  +  p)V-ntp_ 

(eqn  2.14) 


C.  NUMERICAL  SOLUTION 

The  finite-difference  scheme  employed  in  Sankar's  problem  solver  is  based  on  the 
Beam-Warming  approximate  factorization  algorithm  [Ref.  10:  pp.  489-496]  as 
implemented  by  Steger  [Ref.  9].  A  theoretical  description  is  given  in  References  6  and 
7.  Solving  the  governing  equation  in  the  transformed  plane  (Eqn.  2.11)  requires 
determination  of  the  metrics  and  Jacobian,  relating  the  variables  to  corresponding 
quantities  in  the  physical  plane  through  Equations  2.12  and  2.14.  This  is  accomplished 
through  Equations  2.8  and  2.9,  where  central  difference  approximations  are  used  to 
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compute  the  spatial  derivatives  x^,  y^,  x^,  and  at  each  interior  grid  point,  and  xT 
and  yt  are  the  grid  velocity  components.  One-sided  differences  are  used  at  the 
boundaries.  The  governing  equation  (Eqn.  2.11)  is  parabolic  (hyperbolic  when  the 
viscous  terms  are  neglected)  [Ref.  10:  p.  139],  lending  itself  naturally  to  an  implicit 
solution  approach  [Ref  14].  This  allows  the  use  of  a  larger  time  step  than  explicit 
schemes,  and  the  requirements  on  step  size  are  generally  set  by  accuracy  rather  than 
stability  requirements  [Refs.  2,9].  The  use  of  explicit  boundary  conditions  in  the 
present  implementation  introduces  a  step-size  stability  requirement.  Given  an  assumed 
flow  field  at  time  tn,  an  alternating-direction-implicit  (ADI)  procedure  is  used  to 
advance  the  solution  to  time  level  tn  +  J.  Elements  of  the  vector  q*  are  then  given  in 
terms  of  values  at  time  tn  and  increments,  ie.  [q*}n+1  =  {q*}n  +  {Aq*}n  +  1.  In 
Sankar's  program,  the  approximation  method  used  is  an  Euler  implicit  finite-difference 
scheme  [Ref  10:  p.  98],  with  central  spatial  differences  and  backward  time  differences 
for  the  derivatives  in  the  governing  equation.  The  method  is  thus  first-order  time 
accurate  and  second-order  accurate  in  space.  [Ref.  7:  pp.  21-23] 

1.  Linearization  and  Factorization 

Inspection  of  the  terms  of  the  flux  and  viscous  stress  vectors  of  the  governing 
equation  shows  them  to  be  very  non-linear  functions  of  the  unknown  vector  q*.  The 
solution  procedure  involves  the  use  of  truncated  Taylor  series  expansions  to  linearize 
the  equation  [Ref.  10:  pp.  490-491].  The  viscous  terms,  however,  are  treated  explicitly, 
evaluated  from  the  solution  at  the  previous  time  level,  rather  than  following  Steger's 
fully  implicit  treatment  [Ref.  9].  This  leads  to  improved  computational  efficiency  but 
requires  implicit  smoothing  for  stability  at  moderate  to  high  Reynolds  numbers 
[Ref  6].  The  Euler  implicit  form  may  be  obtained  from  Taylor  series  expansion  at  time 
Replacing  the  partial  time  derivative  with  the  backward  difference  operator 

yields: 

[Aq*}n  +  1  *  (q*}n+1  -  {q*}n  =  At(Vt{q*})n  +  1  4-  O(At)2  (eqn  2.15) 

Substituting  Equation  2.11  into  Equation  2.15,  after  replacing  the  spatial  derivatives 
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with  the  central  difference  operators  fij:  and  6^  and  the  time  derivative  with  the 
backward  difference  operator,  produces  the  governing  equation  in  non-linear,  time- 
differenced  form  [Refs.  7,10:  pp.  22-23,490-491]: 

{Aq*}n  +  1  «  -At(6^(F*}n'1'1  +  5n{G*}n+1) 

+  At(S^(R*}n  +  1  +  6n(S*}n+1)  +  O(At)2  (eqn  2.16) 

Taylor  series  expansions  of  the  flux  vectors  are  now  introduced: 

(F*}n  +  1  =  (F*}n  +  [5qF*]n{Aq*}n"t" 1  +  O(At)2 

(G*}n  +  1  .  +  [^G*]n{Aq*}n+1  +  O(At)2  (eqn  2.17) 

Here  [dqF*]n  and  [dqG*]n  are  Jacobian  matrices,  evaluated  under  the  assumption  of  a 
perfect  gas  [Refs.  7,10:  pp.  87-88,491],  which  will  be  denoted  [A]  and  [B]  respectively. 
Substituting  Equation  2.17  into  2.16  and  placing  the  implicit  terms  on  the  left  side  and 
explicit  terms  on  the  right  side  gives  the  linearized  system 

([I]  +  Atd^[A]  +  Atdn[B])  [Aq*}n+1  -  {R}n  (eqn  2.18) 

where 

(R}n  -  -At(6^{F*}n  +  8n(G*}n)  4-  At(5jj{R*}n  +  1  +  8n(S*}n  +  1) 

Since,  as  mentioned  previously,  the  viscous  terms  are  treated  explicitly,  they  have  been 
grouped  together  in  [R}n  with  the  other  terms  to  be  evaluated  at  the  known  time  level 
tn.  [Refs.  6,7:  pp.  23-24] 

Although  linearized,  the  system  of  equations  is  now  in  block  pentadiagonal 
form  and  still  computationally  expensive  to  solve.  Approximate  factorization  is  used  to 
wTite  the  implicit  scheme  as  a  product  of  one-dimensional  factors  so  that  the  problem 
may  be  solved  by  a  sequence  of  relatively  simple  operations,  subject  to  the  additional 
error  introduced  by  the  factorization  approximation.  [Ref.  14]  Factoring  the  left  side 
of  Equation  2.18  then  gives  [Ref.  7:  p.  25] 

([I]  +  ArityA])  ([I]  +  At«n[B])  {Aq*}n  +  1  -  {R}“  (eqn  2.19) 

The  Beam-Warming  algorithm  is  now  implemented  in  three  steps  [Ref.  10:  p.  494]: 


28 


Step  1 

([I]  +  At8^[A])  (Aq**}  -  [R}n 

(eqn  2.20) 

Step  2 

([I]  +  At8n[B])[Aq*}n  +  1  -  {Aq**} 

(eqn  2.21) 

Step  3 

(q*)n  +  1  *  {q*}n  +  {Aq*}n  +  I 

(eqn  2.22) 

Notice  that  steps  one  and  two  represent  sweeps  in  alternating  and  q-  directions, 
leading  to  the  name  alternating-direction-implicit  for  this  method.  Since  the  difference 
operators  5e  and  8^  represent  central  difference  approximations,  they  lead  to  systems 

of  block  tridiagonal  matrix  equations  composed  of  4  *  4  submatrices,  which  are  easily 
solved. 

2.  Application  of  Boundary  Conditions 

On  the  solid  surface  the  flow  tangency  condition  applies  for  the  Euler 
equations,  and  in  addition,  the  no-slip  condition  applies  for  the  viscous  flows 
considered  here.  This  means  the  fluid  and  solid  surface  must  have  the  same  velocity  at 
the  common  boundary,  which  is  to  say  the  contravariant  velocity  components  U  and  V 
(Eqn.  2.13)  are  both  zero.  The  surface  is  assumed  adiabatic  so  dne  -  0.  It  is  also 
assumed,  as  reasonable  approximations  for  high  Reynolds  number  flows,  that  dnp  -  0 
and  5np  =  0  at  the  surface.  Pressure  and  density  at  the  surface  are  then  numerically 
determined  by  two-point  extrapolation,  p.j  «  (4p.  2  -  p.  3)/3  and  p.  j  -  (4p.  - 
pi  3)/3  internal  energy  is  calculated  using  the  relation  p  -  (y  -  lXe  -  0.5p(u2  +  v2)). 
Boundary  conditions  are  specified  following  each  ADI  sweep  for  the  incremental 
quantities  {Aq*}  and  {Aq**}  by  setting  them  equal  to  zero  on  the  solid  surface. 
[Refs.  6,7,11]  The  far-field  boundary  conditions  are  undisturbed  freestream  conditions. 
Conditions  along  the  cut  in  the  C-grid  are  found  by  averaging  the  solution  at  the 
nearest  interior  points  on  either  side  of  the  cut.  [Ref.  7:  pp.  27-29]  The  effect  of  this 
scheme  is  to  update  conditions  at  the  boundaries  explicitly.  The  explicit  boundary 
conditions  are  used  for  their  simplicity  in  the  code,  although  implicit  conditions  would 
generally  be  preferable  for  stability  and  accuracy  of  the  numerical  solution  [Refs.  9,15], 
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3.  Artificial  Dissipation 

When  Taylor  series  expansions  are  substituted  into  the  finite-difference 
approximation  of  a  partial  differential  equation,  and  the  terms  are  rearranged  to 
produce  the  original  differential  equation  plus  a  truncation  error,  the  lowest  order 
truncation  error  term  may  contain  a  second  derivative  that  makes  the  term  similar  in 
appearance  to  the  viscous  terms  in  flow  equations.  Thus  the  use  of  finite-difference 
approximations  introduces  an  "artificial  viscosity"  into  the  solution.  [Ref.  10:  pp. 
89-92]  One-sided  difference  schemes  (upwind  differencing)  are  commonly  used  to 
produce  this  implicit  dissipation.  Additionally,  terms  may  be  added  explicitly  to 
suppress  oscillatory  behavior  of  the  numerical  solution.  Central  difference 
approximations,  given  by  (ui+1  -  u..j)/2Ax,  tend  to  decouple  the  odd  from  the  even 
terms,  since  only  every  other  term  is  used  in  the  difference  formula.  For  example,  the 
sequence  of  nodal  values  ( - 1,+  1,  -  1,+  1,...}  would  give  a  central  difference 
approximation  to  the  first  derivative  of  zero.  [Ref.  14]  Solutions  of  the  Navier-Stokes 
equations  may  "blow  up"  due  to  numerical  oscillations  if  the  mesh  is  not  fine  enough  in 
areas  of  large  pressure  gradients.  To  suppress  high  frequency  numerical  oscillations, 
fourth-order  explicit  dissipation  terms  are  commonly  added  to  algorithms.  [Ref.  10:  pp. 
105,  486]  Physical  dissipation  diffuses  energy,  and  artificial  dissipation  reduces  gradients 
in  the  flow-field  solution  whether  such  diffusion  is  physically  correct  or  numerically 
induced.  The  explicit  method  of  adding  dissipation  has  the  advantage  over  one-sided 
differencing  of  making  the  physical  approximations  clearer  and  permitting  some  control 
over  the  amount  of  non-physical  (numerical)  dissipation.  [Ref.  4]  Although  viscous 
terms  provide  a  dissipative  mechanism,  some  additional  numerical  dissipation  is  ahvays 
necessary;  the  degree  should  be  kept  to  the  minimum  required.  [Refs.  14,10:  p.  92] 

The  Navier-Stokes  code  compiled  by  Sankar  follows  Steger  in  the  use  of 
artificial  dissipation,  with  changes  to  prevent  overshooting  in  the  vicinity  of  shocks. 
[Ref.  16]  Dissipation  is  employed  for  four  main  purposes,  the  first  three  of  which  have 
been  mentioned  previously: 

1.  Suppress  high  frequency  oscillations  in  the  solution  caused  by  the  use  of  central 
differences. 

2.  Correct  for  the  incorrect  initial  conditions,  especially  when  solving  inviscid 
flow's. 

3.  Allow  explicit  treatment  of  viscous  terms. 

4.  Alleviate  restrictive  stability  bounds  on  the  use  of  explicit  damping.  [Ref.  9] 
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The  first  two  reasons  pertain  to  the  right-hand  (explicit)  side  of  Equations  2.19  and 
2.20  and  the  second  two  reasons  to  the  left-hand  (implicit)  side  of  Equations  2.19,  2.20, 
and  2.21.  The  explicit  dissipation  used  is  a  combination  of  second-  and  fourth-order 
terms.  Fourth-order  smoothing  is  normally  used  for  accuracy,  since  second-order 
terms  tend  to  smear  rapid  pressure  variations  near  the  leading  edge.  In  large  pressure 
gradients  such  as  shocks,  however,  these  terms  lead  to  overshooting  in  the  pressure 
distribution.  To  prevent  this,  the  second  derivative  of  pressure  is  used  to  control  the 
degree  of  dissipation,  turning  on  the  second-order  term  and  suppressing  the  fourth- 
order  term  in  the  vicinity  of  shocks  [Ref.  6). 

Upwind  difference  schemes  may  be  rewritten  as  the  sum  of  central  difference 
approximations  to  the  first  and  second  derivatives  as  follows: 

(ui  ~  u..j)/Ax  -  (u.  +  ,  -  u..1)/2Ax  —  (u.  +  1  +  2U-  —  u..,)/2Ax  (eqn  2.23) 

The  second  term  on  the  right-hand  side  may  then  be  regarded  as  an  implicit  dissipation 
term,  giving  a  central  difference  scheme  the  artificial  viscosity  characteristics  of  upwind 
diffeuncing.  Second-order  implicit  smoothing  terms  are  added  to  the  left-hand  side  of 
Equations  2.19,  2.20,  and  2.21,  both  to  allow  the  use  of  explicitly  evaluated  viscous 
terms  and  to  permit  the  use  of  larger  magnitude  explicit  damping  terms  without  loss  of 
stability.  The  dissipation  terms  are  of  at  most  the  same  order  as  the  truncation  errors 
of  the  finite-difference  approximation  involved  [Refs.  6,7:  p.  30-33].  The  coefficients  of 
both  implicit  and  explicit  dissipation  terms,  Ej  and  ee,  are  user  selectable  options  that 
control  the  magnitude  of  the  artificial  viscosity. 

D.  TURBULENCE  MODEL 

The  Reynolds  form  of  the  Navier-Stokes  equations  introduces  new  unknowns  in 
the  form  of  turbulent  stress  and  heat-flux  terms,  requiring  additional  relations  to  make 
solution  of  the  system  possible.  This  is  known  as  the  closure  problem,  and  the 
approach  usually  taken  to  resolve  it  is  turbulence  modeling.  [Ref.  10:  pp.  207-208, 
221-235]  Under  the  Boussinesq  assumption,  the  coefficients  of  viscosity  and  thermal 
conductivity  are  replaced  by  the  combinations  |i  +  and  k  +  kt  representing  the  sum 
of  laminar  and  turbulent  components.  In  practice  the  thermal  conductivity  is  usually 
determined  from  the  relationships  k  =  cpp/Pr  and  kt  -  cppt/Prt.  In  the  solver  used 
here  a  constant  value  of  Pr  =  1  is  assumed,  and  total  viscosity  is  used  in  Equation  2.3 
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rather  than  directly  calculating  the  conductivity.  For  the  present  calculations,  the 
viscous  work  and  conductivity  terms  R4  and  S4  (Eqn.  2.5)  were  set  equal  to  zero.  The 
laminar  viscosity  is  easily  found  from  the  Reynolds  number  and  freestream  velocity, 
and  an  algebraic  turbulence  model  is  used  to  determine  the  turbulent  eddy  viscosity. 

Turbulence  modeling  lies  at  the  frontier  of  research  in  computational  fluid 
dynamics.  The  use  of  Navier-Stokes  solvers  as  practical,  predictive  tools  depends  on 
the  development  of  adequate  turbulence  models.  For  many  uses,  however,  the 
dissipative,  diffusive  nature  of  turbulence  has  enabled  simple  algebraic  models  to 
perform  beyond  their  theoretical  limitations  [Ref.  17].  In  the  present  case  the  complex, 
history-dependent  viscous  effects  dominating  the  flow  are  clearly  beyond  the  capacity 
of  simple  algebraic  models  based  strictly  on  local  flow  properties.  Unfortunately,  there 
are  no  completely  satisfactory  alternatives.  The  Baldwin-Lomax  model  was  chosen  for 
simplicity,  while  acknowledging  that  it  might  be  unsuitable  for  the  massively  separated 
flows  experienced  in  dynamic  stall  [Ref.  6].  When  it  was  discovered  that  the  maximum 
lift  was  underpredicted,  a  modification  was  incorporated  which  will  be  discussed 
following  a  description  of  the  original  Baldwin-Lomax  model. 

1.  Baldwin-Lomax  Model 

The  Baldwin-Lomax  model  is  based  on  Cebeci's  two-layer  model,  with  pt 
given  by  the  inner-layer  formulation  out  to  the  smallest  normal  distance  at  which  the 
inner-  and  outer-layer  formulas  give  the  same  value  for  the  viscosity.  In  the  inner  layer 
a  simple  mixing  length  formula  is  used  with  the  eddy  viscosity  proportional  to  the 
magnitude  of  the  local  vorticity  times  the  square  of  the  normal  distance  from  the 
surface: 

(i>tW  =  Pf2l“l 

w’here  the  mixing  length 

£  =  KyD  (eqn  2.24) 


and  D  is  the  Van  Driest  damping  function,  given  by: 
D  =  [1  -exp(-y  +  /A  +  )] 


(eqn  2.25) 
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The  damping  constant  A+  ■  26,  and  the  von  Karman  constant,  k,  is  taken  as  0.4.  In 
the  outer  layer,  the  locally  constant  eddy  viscosity  is  given  by: 


^T^outer  "  ^cpP^wake^Kleb^) 


(eqn  2.26) 


The  key  feature  is  the  definition  of  a  function 

Fwake  =  1™n(ymaxFmax’  ^k^max*-1  dif^Fmax)  (ecln  2.27) 

where  Fmax  is  the  maximum  of  the  function 

F(y)  -  yl«|D  (eqn  2.28) 

arM*  Tmax  *s  va*ue  of  y  at  which  Fmax  occurs.  The  KlebanofF  intermitancy 
function  FKleb,  has  the  effect  of  causing  the  eddy  viscosity  to  approach  zero  in  the  far 
field.  Udif  is  the  difference  between  the  maximum  and  minimum  velocities  in  the 
velocity  profile.  The  values  of  the  constants  used  are  K  -  0.0168,  C£  ■  1.6,  and 
Cwk  -  0.3.  [Refs.  2,6,7:  pp.  16-19] 

2.  Modified  Model 

It  has  been  observed  that  algebraic  turbulence  models  of  the  Cebeci-Smith 
type  [Ref.  10:  p.  225]  perform  well  in  attached  flows  at  moderate  angles  of  attack  but 
overpredict  both  the  rise  in  turbulent  shear  stress  in  adverse  pressure  gradients,  and  the 
stress  decrease  when  pressure  gradients  are  relieved  [Ref.  18].  The  Baldwin- Lomax 
model  was  found  by  Sankar  to  cause  an  early  prediction  of  flow*  separation  at  high 
angles  of  attack.  It  was  assumed  that  this  was  due  to  underprediction  of  the  turbulent 
length  scale  at  high  angles.  The  function  F,  defined  in  Equation  2.28,  has  in  some 
cases  been  found  to  possess  a  number  of  local  maxima,  leading  to  large  length  scale 
variations  depending  on  the  maximum  chosen.  This  presented  a  possibility  for  simple 
modification  to  increase  the  outer-layer  velocity  and  length  scales.  The  definition  of 
Fmax  was  defined  to  be  the  value  of  the  function  F  w’hen  the  product  yF(y)  is 
maximum.  In  preliminary  investigation  this  change  was  found  to  increase  stall  angles 
but  to  lead  to  premature  reattachment  during  the  downstroke  in  dynamic  calculations. 
For  this  reason  its  use  is  recommended  only  until  stall  onset.  [Ref.  19] 
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III.  DESCRIPTION  OF  THE  CODE 


The  version  of  Sankar's  code  used  in  the  present  study  is  presented  in  Appendix 
A.  Appendix  B  contains  instructions  and  guidelines  for  running  the  program. 
Modifications  were  made  in  the  output  scheme  to  produce  output  of  pressure 
coefficients  and  flow-field  plotting  data  at  a  specified  number  of  equal  time  intervals 
during  each  cycle,  and  job  cards  were  prepared  for  save,  restart,  and  various  output 
options.  Other  features  of  the  program  were  not  altered,  although  comments  were 
added.  Subroutines  in  the  program  supplied  by  Sankar  were  fully  vectorized  for  the 
Cray  computer  with  the  exception  of  subroutine  EDDY.  Execution  times  on  the 
NASA  Ames  Research  Center's  Cray  X-MP/48  computer,  for  full  viscous  dynamic 
calculations,  are  approximately  0.318  seconds  per  time  step. 

All  program  variables  are  nondimensionalized.  This  has  the  effect  of  multiplying 
the  right  side  of  the  governing  Equation  2.11  by  a  factor  of  Re'1  but  does  not  affect 
the  basic  form  of  the  equation  (Ref.  10:  pp.  191-193J.  The  nondimensionalization  is 
effected  as  follows: 

1.  density  with  respect  to  px 

2.  length  with  respect  to  c 

3.  time  with  respect  to  c/a  ^ 

4.  velocities  with  respect  to  a^ 

5.  total  energy  with  respect  to  Pooaoo2 

Under  this  scheme  the  number  of  time  steps  required  to  complete  an  oscillation  cycle  is 
given  by  Jt/(k  x  x  At). 

A.  MAIN  PROGRAM 

The  program  begins  execution  by  reading  the  input  data.  This  consists  of  the 
dynamic  stall  parameters,  <*0,  ctj,  k,  M,  and  Re;  program  control  information,  including 
time-step  size,  explicit  viscosity  coefficient,  and  flags  for  features  such  as  the  modified 
Baldwin-Lomax  turbulence  model;  and  the  airfoil  coordinates  and  related  information. 
Table  1  summarizes  the  subroutine  calls.  The  grid  generation  routine,  AIRFOL,  is 
called,  and  for  viscous  flows,  CLUSTR  is  then  called  to  recompute  the  constant-ij  grid 
lines  so  that  the  first  point  from  the  airfoil  surface  is  at  the  user-specified  distance.  The 
initial  conditions  are  established  next,  and  if  the  program  is  being  restarted  from 
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previous  calculations,  the  values  of  time  and  the  unknown  flow-field  vector  q*  are 
replaced  with  those  now  read  from  logical  unit  7.  Finally,  subroutine  METRIC  is 
called  to  compute  the  metrics  and  transformation  Jacobian,  and  the  iterative 
calculations  begin. 


TABLE  1 

SUBROUTINE  CALLS  FROM  THE  MAIN  PROGRAM 

Program  MAIN 

AIRFOL 

-  generate  grid 

CLUSTR 

-  cluster  grid  points  near  the  surface  for  viscous  calculations 

METRIC 

-  compute  metrics  and  Jacobian 

Loop: 

ROTGRID 

-  rotate  grid  to  current  angle  of  attack 

METRIC 

-  recompute  metrics 

SLPS 

-  perform  ADI  sweeps 

WALLBC 

-  apply  boundary  conditions  on  surface  and  cut 

CPPLOT 

-  output  surface  pressure  distribution 

LOAD 

-  compute  integrated  lift,  moment,  and  drag  coefficients 

The  iterative  solution  loop  is  designed  to  execute  the  ADI  sweeps  for  a  number 
of  steps  specified  in  the  input  data.  For  restarts  of  dynamic  calculations,  the  iterative 
solution  loop  begins  its  first  cycle  by  rotating  the  grid  (subroutine  ROTGRID)  based 
on  the  time  read  from  the  previously  saved  solution  and  the  frequency  of  oscillation. 
No  rotation  is  required  at  the  initiation  of  dynamic  calculations  or  for  steady-state 
calculations  since  the  initial  conditions  set  the  freestream  flow  direction  at  an  angle 
equal  to  the  initial  angle  of  attack.  On  each  subsequent  iteration  the  grid  is  rotated  an 
incremental  amount  based  on  the  size  of  the  time  step.  Following  each  grid  rotation, 
subroutine  METRIC  is  again  called  to  recompute  the  metrics  since  the  computational 
plane  is  stationary.  The  actual  solution  is  now  accomplished  by  a  call  to  SLPS,  and 
boundary  conditions  are  applied  at  each  step  by  WALLBC.  The  completed  solution 
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may  now  be  used  to  produce  normal  program  output,  at  intervals  specified  in  the 
program.  CPPLOT  is  called  to  output  the  surface  pressure  distribution,  and  LOAD 
computes  the  integrated  lift,  drag,  and  moment  coefficients.  After  exiting  the  loop, 
various  output  options  may  be  selected,  including  the  option  of  saving  the  current 
solution  for  a  subsequent  restart;  creating  the  plotting  files;  and  writing  the  complete 
flow-field  solution,  including  the  velocity  magnitude,  eddy  viscosity,  and  normal 
distance  from  the  surface  at  each  grid  point.  The  present  version  of  the  program  is 
modified  to  exit  at  a  number  of  equal  time  intervals  each  cycle,  as  specified  within  the 
program.  This  was  done  to  generate  data  files  for  PLOT3D  (an  Ames  Research  Center 
plotting  code)  at  equal  phase  intervals,  and  also  to  save  the  current  solution  so  it  may 
be  written  to  Cray  tapes  for  possible  future  use. 

B.  GRID  GENERATION 

The  call  to  subroutine  AIRFOL  initiates  the  basic  grid  generation  process. 
AIRFOL  begins  by  reading  the  airfoil  shape  input  parameters.  These  include  the 
number  of  points  on  the  upper  and  lower  surfaces  and  a  flag  indicating  airfoil 
symmetry.  For  symmetric  airfoils,  only  the  upper  surface  coordinates  are  read,  and  the 
lower  surface  coordinates  are  defined  as  (XL,  YL)  =  (XU,  -YU).  Then  the 
coordinates  are  scaled  with  respect  to  the  airfoil  chord  length,  c,  and  the  actual  grid 
generation  begins.  The  sequence  of  subroutine  calls  is  contained  in  Table  2. 

The  first  step  in  the  grid  generation  is  to  fill  in  the  definition  of  the  airfoil 
surface.  Subroutine  SING  is  called  to  determine  the  singular  point  (x^-,  y^;),  defined 
earlier  as  lying  midway  between  the  nose  and  the  origin  of  the  leading-edge  radius. 
SING  applies  a  three-point  parabolic  curve  fit  to  the  leading-edge  points  to  determine 
the  radius  of  curvature.  The  angle  of  the  trailing-edge  upper  surface  relative  to  the 
leading  edge  and  the  average  of  upper-surface  leading-  and  trailing-edge  slopes  are  also 
computed  for  use  in  determining  the  wake  shape.  AIRFOL  next  calls  subroutine 
TABINT  to  compute  the  additional  points  required.  TABINT  first  performs  the 
mapping  to  an  intermediate  plane  described  by  Equation  2.1.  The  distance  required 
between  points  to  place  97  points  on  the  airfoil  surface  is  calculated.  Then  the  actual 
determination  of  the  points  is  performed  by  a  standard  polynomial  interpolation 
routine,  TAINT.  Finally,  TABINT  maps  the  smoothed  surface  back  to  the  physical 
plane. 
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TABLE  2 

GRID  GENERATION  SUBROUTINES  CALLED 

AIRFOL 

-  called  by  MAIN 

SING 

-  determine  singular  point  (x^r,  y^) 

TABINT 

-  find  additional  points  on  airfoil  surface 

TAINT  -  interpolation  routine 

WRAP 

-  calculate  stretching  in  *}-direction 

CLUSTR 

-  called  by  MAIN  for  viscous  flows 

STRTCH 

-  calculate  new  stretching  factor 

TAINT 

-  locate  new  normal  grid  point  locations 

A  smooth  airfoil  shape  has  now  been  defined.  It  remains  only  to  determine  the 
shape  of  the  wake— the  'cut"  in  the  physical  domain— before  constructing  the  grid. 
AIRFOL  uses  the  trailing-edge  angle  computed  by  SING  to  calculate  a  wake  shape 
allowing  the  flow  to  leave  the  trailing  edge  smoothly.  This  cut  generally  corresponds 
to  the  tangent  of  the  mean  camber  line,  and  for  symmetric  airfoils,  it  is  just  the 
extended  chord  line.  [Ref.  7:  p.  41]  It  remains  at  a  fixed  angle,  and  the  entire  grid  is 
rotated  along  with  the  airfoil.  (See  Fig.  2.4) 

Construction  of  the  grid  is  finally  completed  by  subroutine  WRAP.  The 
stretching  in  the  normal  direction  is  computed.  Then  the  airfoil  and  wake  are 
unwrapped  using  Equation  2.1.  This  creates  a  grid  consisting  of  constant-*}  lines 
equidistant  from  the  unwrapped  surface  and  constant-^  lines  parallel  to  the  »}-axis. 
This  grid  is  used  for  inviscid  (Euler)  calculations.  Following  the  return  from  WRAP, 
the  last  step  performed  by  AIRFOL  is  to  map  the  entire  grid  back  into  the  x-y  plane, 
and  the  coordinate  axis  is  shifted  to  the  quarter-chord  point. 

If  the  Reynolds  number  is  greater  than  zero,  program  MAIN  calls  subroutine 
CLUSTR  to  construct  new  constant-*}  lines.  (Setting  the  Reynolds  number  to  zero  or 
a  negative  value  is  a  flag  for  inviscid  calculations).  For  each  grid  x-coordinate,  the  old 
stretching  in  the  normal  direction  is  computed;  that  is,  the  physical  distance  from  the 
surface  to  each  point  along  constant-^  lines  is  calculated.  Then  subroutine  STRTCH 
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is  called  to  calculate  the  new  stretching.  In  STRTCH  the  user-input  distance  of  the 
first  point  off  the  wall  and  the  distance  just  calculated  by  CLUSTR  to  the  farthest  grid 
point  are  used  as  the  inner  and  outer  bounds  of  the  grid.  A  stretching  factor  R  is 
desired  to  give  a  geometric  progression  of  grid  spacings  such  that,  given  the  distance  of 
the  first  point  off  the  wall,  the  ratio  between  successive  spacings  is  R.  This  correct 
stretching  factor  is  calculated  by  iteration  using  Newton's  method.  Once  the  stretching 
factor  has  been  determined,  CLUSTR  calls  the  interpolation  routine  TAINT  to  locate 
the  new  grid  coordinates. 

Construction  of  the  grid  in  both  the  physical  and  computational  planes  is  now 
complete.  The  physical  grid  is  easily  rotated  by  subroutine  ROTGRID  using  the 
relations  x  *  x  cos  (Aa)  —  y  sin  (Aa)  and  y  =  y  cos  (Ao)  +  x  sin  (Aa).  The  only 
step  remaining  is  the  calculation  of  the  transformation  metrics  and  Jacobian. 
Subroutine  METRIC  is  called  by  MAIN  after  construction  of  the  grid  and  after  each 
rotation.  METRIC  performs  calculations  as  outlined  in  the  previous  chapter  under 
"Numerical  Solutions".  The  grid  velocity  components  are  xT  and  yr  determined  from 
the  angular  velocity  and  grid  coordinates  in  the  physical  plane  (since  coordinates  were 
defined  with  respect  to  the  quarter  chord).  Central  differences  are  used  at  interior 
points  and  one-sided  differences  at  the  boundaries-two-point  downstream  and  three- 
point  elsewhere-to  compute  x^,  x^,  y^,  and  y^.  This  calculation  is  simplified  by  the 
assumption  of  a  unit  grid  spacing  in  the  computational  plane.  Finally,  the  relations 
given  by  Equations  2.8  and  2.9  are  used  to  compute  the  Jacobian  and  metrics. 

C.  COMPUTATIONAL  STEPS 

The  basic  features  of  the  computational  scheme  were  outlined  in  the  preceding 
chapter.  The  organization  of  subroutines  to  implement  this  scheme  is  presented  in 
Table  3.  Program  MAIN  calls  subroutine  SLPS  to  perform  the  primary  calculation 
steps.  SLPS  controls  the  execution  of  the  ADI  sweeps,  assembles  the  matrices,  collects 
the  damping  terms,  and  calls  the  matrix  solvers.  The  value  of  the  implicit  damping 
coefficient  is  set  in  relation  to  the  explicit  coefficient.  In  the  present  implementation  £j 
=  20eg.  If  the  reduced  frequency  is  less  than  .001,  a  local  time-step  option  is  used  to 
accelerate  convergence  to  a  steady-state  solution  by  basing  the  time  step  on  the  local 
flow  conditions.  The  calculations  begin  by  calling  subroutine  DISSIP  to  compute  the 
explicit  dissipation  terms  described  in  the  last  chapter  and  add  them  to  the  right  side  of 
the  governing  equation  (Eqn.  2.20).  In  the  ^-direction  DISSIP  uses  a  blend  of  second- 
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and  fourth-order  terms.  A  switching  function  is  computed,  based  on  the  second 
derivative  of  pressure.  This  function  causes  the  fourth-order  terms  to  be  set  to  zero 
when  the  pressure  gradient  exceeds  a  value  specified  within  the  program.  Otherwise 
the  second-order  terms  are  set  to  zero.  The  large  pressure  gradients  experienced  in  the 
vicinity  of  shocks  are  not  expected  in  the  q-direction,  so  only  fourth-order  terms  are 
added.  The  arrays  DQ1-DQ4  contain  the  combined  dissipation  terms. 


TABLE  3 

SUBROUTINES  USED  IN  COMPUTATIONAL  STEPS 

SLPS 

-  called  by  MAIN 

DISSIP 

-  compute  explicit  dissipation  terms 

RESI 

-  compute  inviscid  right-hand  (explicit)  terms 

STRESS 

-  compute  viscous  right-hand  terms 

EDDY  -  determine  viscosity  coefficient 

AMAT1 

-  compute  Jacobian  matrix  [A] 

MATRX1 

-  invert  assembled  matrix  in  the  ^-direction 

AMAT2 

-  compute  Jacobian  matrix  [B] 

MATRX2 

-  invert  assembled  matrix  in  the  q-direction 

WALLBC 

-  called  by  MAIN  to  apply  surface  boundary  conditions 

Subroutine  SLPS  next  calls  RESI  to  compute  the  inviscid  right-hand  terms  at  the 
known  time  level  (Eqn.  2.18).  Working  first  with  the  ^-derivative  flux  terms,  RESI 
begins  by  calculating  the  contravariant  velocity  component  U,  in  the  computational 
plane  (Eqn.  2.13).  The  vector  F*  (Eqn.  2.14)  is  then  formed  using  the  contravariant 
velocity.  Then  -At[^(F*)n]  is  computed  using  standard  central  differences  and  added 
to  the  arrays  DQ1-DQ4.  These  steps  are  repeated  in  the  q-direction,  using  the 
contravariant  velocity  component  V,  and  -At[dn(G*)n]  is  added  to  the  DQ  arrays. 

For  viscous  flows,  subroutine  STRESS  is  called  next  to  compute  the  viscous 
terms  and  add  them  to  the  right-hand  side  of  the  equation.  The  first  step  in  STRESS 
is  a  call  to  EDDY  to  compute  the  total  (turbulent  and  laminar)  viscosity.  (EDDY  is 
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not  called  on  the  first  ten  steps  of  the  initial  dynamic  run  to  allow  the  solution  to  settle 
a  bit).  The  eddy  viscosity  routine  begins  by  initializing  the  viscosity  throughout  the 
field  to  the  laminar  value  calculated  from  the  Mach  number  and  Reynolds  number.  If 
laminar  calculations  are  desired,  the  subroutine  should  be  exited  at  this  point.  A 
completely  turbulent  flow  field  is  assumed,  so  no  transition  point  must  be  determined. 
The  subroutine  works  with  one  grid  ^-location  at  a  time.  The  calculations  performed 
were  described  in  the  last  chapter  under  "Turbulence  Model".  The  variables  defined  in 
that  section  are  initialized  to  low  values  which  will  be  reset  as  the  flow  field  is  scanned 
in  a  normal  direction  from  the  surface.  Beginning  at  the  airfoil  boundary,  the 
derivatives  u^,  v^  are  calculated  by  one-sided  differences,  and  the  metrics  and  Jacobian 
are  again  computed  as  described  previously.  These  values  are  all  local  to  subroutine 
EDDY.  The  no-slip  condition  is  now  applied  to  determine  the  viscous  stress,  tw,  and 
skin  friction.  The  friction  value,  CF,  is  one  of  the  selectable  normal  output  values 
from  MAIN.  The  Van  Driest  damping  factor,  y  +  /A  +  ,  is  calculated  and  the 
subroutine  begins  scanning  the  grid  points  outward  from  the  surface.  The  values  uc, 
v^,  u^,  v^:  the  metrics  and  Jacobian;  and  y,  and  F  are  calculated,  and  the  values 
of  Fmax  and  ymax  are  found  (Eqns.  2.27  and  2.28).  If  the  angle  of  attack  is  less  than 
the  user-specified  value  ALFAI  and  is  increasing  or  constant,  the  modified  turbulence 
model  is  used  and  F—.,,  is  the  value  of  F  where  F  *  y  is  maximum.  Above  ALFAI  on 
the  upstroke  and  on  the  entire  downstroke,  the  original  Baldwin-Lomax  model  is  used 
and  Fmax  is  the  maximum  value  of  F.  Then  the  length  scale,  C,  is  found  using  the  Van 
Driest  function,  and  the  inner-layer  viscosity  is  computed  (Eqn.  2.24).  Next,  the  outer- 
layer  formula  previously  described  is  used  to  find  a  value  for  viscosity  at  each  normal 
grid  location;  and  the  point  where  the  inner-  and  outer-layer  values  first  cross  is  found. 
The  final  step  performed  by  EDDY  is  to  add  the  laminar  viscosity  to  the  appropriate 
inner-  or  outer-layer  turbulent  viscosity. 

Now  that  the  viscosity  has  been  calculated,  subroutine  STRESS  can  carry  out 
evaluation  of  (d^R*  +  d^S*)  at  the  known  time  level,  where  R*  and  S*  are,  as  given 
by  Equation  2.12,  combinations  of  the  vector  R  and  S.  First  the  derivatives  in  the  In¬ 
direction  are  computed  by  central  differences.  A  more  efficient  discretization  of  the 
spatial  derivatives  of  viscous  terms  is  realized  by  evaluation  at  half  points,  rather  than 
at  the  nodal  points.  This  reduces  the  number  of  nodal  points  involved  from  five  to 
three.  [Ref.  7:  p.  26]  Therefore,  STRESS  computes  derivatives,  metrics,  Jacobian,  and 
an  average  viscosity  at  the  half  points.  These  values  are  used  to  determine  txx,  rxy, 
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v  R4’  and  S4-  In  the  present  implementation  R4  and  S4,  the  viscous  work  (energy 
dissipation)  and  conduction  (diffusion)  terms,  are  set  equal  to  zero.  Finally  the 
differenced  value  of  d ^  R*  is  added  to  the  DQ1-DQ4  arrays.  The  same  steps  are  then 
followed  for  the  ^-derivative  terms. 

Subroutine  SLPS  has  now  assembled  in  the  arrays  DQ1-DQ4,  the  elements  of  the 
right-hand  side  of  Equation  2.20  given  by  Equation  2.18,  with  explicit  damping  added 
but  without  the  factor  At.  Next,  the  subroutine  AM  AT  l  is  called  to  compute  the 
Jacobian  matrix  [A]  =  d^F*.  and  the  left-hand  side  of  Equation  2.20  is  then  assembled. 
The  second-order  implicit  damping  terms  are  computed  and  added  to  the  left  side. 
Finally,  the  right-hand  side  (DQ  arrays)  is  multiplied  by  At  and  stored  in  matrix  GG, 
giving  the  final  form  of  the  right-hand  side  of  Equation  2.20.  A  standard  matrix 
inversion  routine,  MATRX1,  obtains  the  solution  to  Equation  2.20,  and  Step  1  of  the 
Beam- Warming  algorithm,  the  ^-sweep,  is  now  complete.  The  solution,  corresponding 
to  Aq**  in  Equation  2.20,  is  now  stored  in  the  DQ1-DQ4  arrays  to  become  the  right- 
hand  side  of  Equation  2.21.  The  steps  used  in  the  ^-sweep  are  now  repeated  for  the  r|- 
direction.  AMAT2  computes  the  Jacobian  matrix  [B]  =  5qG*,  and  the  left-hand  side 
of  Equation  2.21  is  assembled  with  implicit  damping  again  added.  MATRX2  is  called 
to  obtain  the  solution  {Aq*}11  ^  and  complete  Step  2  of  the  Beam- Warming  algorithm. 
Step  3  is  accomplished  by  updating  the  flow  variables  (Eqn.  2.22).  Values  at  the 
boundaries  are  not  updated.  The  outer  boundaries  remain  at  undisturbed  freestream 
values  making  a  specific  step  to  apply  boundary  conditions  unnecessary.  This 
completes  the  solution  procedure,  but  as  a  monitoring  feature,  the  density  elements  of 
Aq*  (now’  stored  in  DQ1)  are  scanned  for  the  largest  value.  This  value,  along  with  its 
grid  coordinates,  is  output  at  specified  intervals,  giving  an  indication  of  the  most 
rapidly  changing  area  of  the  flow  field  and  the  magnitude  of  the  density  change  there. 

A  solution  has  now  been  obtained  by  program  MAIN.  It  remains  only  to  apply 
the  surface  boundary  conditions.  This  is  accomplished  by  subroutine  WALLBC. 
Boundary  conditions  were  discussed  in  the  last  chapter,  and  the  implementation  is 
straightforward.  Conditions  along  the  cut  are  uveraged  from  above  and  below.  The 
contravariant  velocity  component  U  at  the  two  points  nearest  to  the  surface  is  used  to 
determine  U  at  the  surface  for  Euler  calculations;  for  viscous  flows  U  -  0  at  the 
surface.  The  appearance  of  the  calculation  steps  is  complicated  slightly  by  use  of  a 
scheme  that  is  applicable  to  both  viscous  and  inviscid  flows.  Surface  flow’  properties 
are  finally  obtained  by  extrapolating  from  interior  points. 
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IV.  RESULTS  AND  DISCUSSION 


Nine  simulations  are  described  in  this  chapter.  Several  variations  were  tried  to 
investigate  the  effects  of  varying  input  parameters,  some  of  which  are  mentioned 
without  detailed  presentation  of  results.  The  intention  was  to  obtain  complete 
solutions  under  a  range  of  conditions  of  dynamic  effects,  Reynolds  number,  Mach 
number,  and  turbulence  model,  whether  or  not  such  solutions  are  the  most  accurate 
obtainable  under  those  conditions.  It  was  hoped  thereby  to  gain  some  understanding 
of  the  most  effective  employment  of  the  code,  its  limitations,  and  possibilities  for 
improving  its  performance.  Table  4  lists  the  input  parameters  for  the  simulations.  The 
calculations  were  grouped  into  the  following  areas: 

1.  Steady-state  calculations  at  and  beyond  maximum  lift. 

2.  Comparison  with  deep  stall  experimental  data,  using  the  original  Baldwin- 
Lomax  turbulence  model  and  modification. 

3.  Comparison  to  experimental  data  for  fully  attached  flow  using  both  turbulence 
model  versions. 

4.  Low  Reynolds  number  predictions  at  Mach  numbers  of  .3  and  .5. 

All  of  the  runs  were  made  with  the  implicit  damping  coefficient  £]  at  a  value  of  twenty 
times  the  explicit  value  £^  (input  as  "WW").  The  distance  of  the  first  point  off  the 
solid  boundary  was  set  to  0.00005  for  all  but  the  steady-state  calculations  beyond 
maximum  lift. 


Three  basic  types  of  output  were  obtained: 

1.  The  flow- field  solution  values  p,  pu,  pv,  and  e  and  the  corresponding  grid 
coordinates  at  twelve  eaual  time  intervals  for  each  cycle,  for  use  in  the  plotting 
program  PLOT3D  bv  Pieter  Buning  of  NASA  Ames.  This  routine  was  used  to 
produce  contour  plots  of  pressure,  density,  Mach  number,  and  stream  function. 

2.  Surface  pressure  coefficients  at  96  eaual  time  intervals  for  each  cycle.  This  data 

wras  used  in  the  plotting  routine  CARPET,  written  by  Rosalie  Lefkowitz  of 
DfssMjv  to^  produce  pressure  distribution  carpet  plots  utilizing  the 

3.  The  normal  text  output  file  containing,  for  96  poipts  each  cycle,  the  integrated 
lft,  drag,  and  moment  coefficients;  pressure  distributions;  the  values  aqd 
ocations  of  the  largest  solution  increments  after  every  ten  time  steps:  and  skin 
riction  values.  These  coefficients  w'ere  plotted  using  the  D1SSPLA  plotting 
package  QPLOT  available  on  the  NASA  Ames  VAX  computers. 

The  use  of  equal  time  intervals  for  output  clusters  the  data  near  maximum  and 

minimum  angles  of  attack,  where  interest  is  often  greatest.  Tht  solution  uras  also 

saved  on  Cray  tapes,  at  twelve  points  each  cycle,  for  future  program  restarts,  so  a  more 
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thorough  investigation  can  be  conducted  of  any  interval  with  minimal  use  of  computer 
time. 


TABLE  4 

INPUT  CONDITIONS  FOR  COMPUTER  SIMULATIONS 

1.  a  =  12°,  Moo  =  -3,  Re  -  1,000,000 

Baldwin-Lomax  Turbulence  Model 

2*  a  ■  14°,  Mqo  -  .3,  Re  -  1,000,000 

Baldwin-Lomax  Turbulence  Model 

3.  a  .«  13.5°,  Moo  *  -301,  Re  -  3,910,000 

Modified  Turbulence  Model 

4.  aQ  -  15°,  -  10°,  Moo  -  -283,  Re  -  3,450,000 

Baldwin-Lomax  Turbulence  Model 

5.  a0  =  15°,  a,  -  10°,  Mqo  -  .283,  Re  -  3,450,000 

Modified  Turbulence  Model 

6.  aQ  -  8°,  cr j  -  10°,  Moo  “  -184,  Re  -  2,450,000 

Baldwin-Lomax  Turbulence  Model 

7.  a0  *  8°,  at  -  10°,  M*  -  .184,  Re  -  2,450,000 

Modified  Turbulence  Model 

8.  a0  -  15°,  aj  -  10*,  Moo  *  -284,  Re  -  345,000 

Baldwin-Lomax  Turbulence  Model 

9.  a0  »  15°,  a,  -  10°,  Mqo  -  -5,  Re  -  345,000 

Baldwin-Lomax  Turbulence  Model 
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A.  STEADY-STATE  CALCULATIONS 

Before  beginning  dynamic  simulations,  the  performance  of  the  program  in 
constant  angle  of  attack  calculations  was  explored.  A  verification  series  of  runs  at 
increasing  angles  of  attack  was  made  using  the  conditions  reported  by  Sankar  in 
Reference  6.  The  results  agree  with  Sankar's  and  are  not  reproduced  here.  The  most 
notable  discrepancy  from  experimental  results  is  the  underprediction  of  maximum  lift 
using  the  Baldwin- Lomax  turbulence  model,  with  separation  occurring  between  one 
and  two  degrees  below  experimental  values.  The  local  time-step  option,  incorporated 
since  Sankar's  results  were  published,  produced  converged  solutions  in  fewer  than  1500 
steps  at  even  the  higher  angles  of  attack.  In  addition  to  the  integrated  coefficients  and 
surface  pressure  distribution  reported  by  Sankar,  the  flow  field  was  graphically 
analyzed.  Figure  4.1  contains  density,  pressure,  Mach  number,  and  stream  function 
plots  at  12  degrees  angle  of  attack. 

Next  the  behavior  of  the  solution  above  static  stall  angles  was  investigated.  The 
original  Baldwin-Lomax  model  was  used  with  a  Mach  number  of  .3  and  a  Reynolds 
number  of  one  million  at  14  degrees  angle  of  attack.  The  time-accurate  mode  was 
invoked  by  using  a  reduced  frequency  of  .002  with  zero  amplitude  of  oscillation.  The 
step  size  was  .005  and  the  distance  to  the  first  point  off  the  wall  was  .0001.  Figure  4.2 
presents  contour  plots  of  density  and  pressure  coefficient  after  5100,  5610,  and  6120 
steps.  Based  on  the  Mach  number  and  time  step,  this  should  represent  a  total 
movement  of  three  chord  lengths  of  a  particle  in  the  freestream.  The  vortex,  however, 
remains  on  the  airfoil  and  gradually  diminishes  in  intensity.  The  integrated  lift 
coefficient  steadily  decreased  throughout  the  run  to  a  value  of  less  than  0.6. 

The  modified  turbulence  model  was  applied  to  the  conditions  of  frame  4019  of 
Reference  8  (volume  3,  page  46)  at  the  experimentally  observed  maximum  lift  angle  of 
attack  of  13.5  degrees.  The  modified  turbulence  model  provides  excellent  results.  For 
a  Reynolds  number  of  3.91  million  and  Mach  number  of  .301,  the  local  time-step 
option  was  used  with  an  explicit  viscosity  input  of  WW  =  2.  After  1500  time  steps  the 
experimental  maximum  lift  coefficient  of  1.4  w’as  obtained,  and  the  peak  pressure 
coefficient  of  8.1  appears  to  be  very  close  to  the  experimental  value.  Moment  and  drag 
coefficient  values  also  seem  reasonable.  The  theoretical  pressure  distribution  at 
maximum  lift  is  shown  in  Figure  4.3,  compared  to  the  experimental  distribution  taken 
from  Reference  8. 
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Figure  4. 1  Steady-State  Attached  Flow 
Mw-.3,  Re-  1.00 x  to6,  a -12° 

Top-Density  (1),  Pressure  (r)  Bottom-Mach  No.  (1),  Stream  Function  (r). 
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Figure  4.2  Steady-State  Separated  Flow 
M x)  =  .3,  Re=  1.00  x  106,  o=  14° 

Density  (1)  and  Pressure  (r)  at  5100,  5610,  and  6120  time  steps. 
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PRESSURE  COEFFICIENT 


Figure  4.3  Steady-State  Pressure  Distribution  at  Maximum  Lift 
Moo -.301,  Re- 3.91  x  I06,a-13.5# 

Inset- Experimental  Data  (Ref.  8). 
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B.  COMPARISON  WITH  DEEP  STALL  EXPERIMENTAL  DATA 

The  first  dynamic  case  studied  was  under  the  experimental  conditions  of  frame 
9218  of  Reference  8  (volume  3,  page  146).  These  are  the  same  conditions  as  Sankar 
used  for  comparison  in  Reference  6.  Sankar  reported  only  the  integrated  lift,  moment, 
and  drag  coefficients,  however,  and  here  the  surface  pressure  distributions  and  the  fiow 
field  about  the  airfoil  were  investigated.  Present  calculations  also  used  the  exact  Mach 
number  and  full  Reynolds  number,  and  runs  were  made  with  both  the  original  Baldwin- 
Lomax  model  and  the  modification.  This  experimental  case  closely  corresponds  to  the 
model  dynamic  stall  (Figure  1.1),  showing  the  formation  of  a  strong  vortex,  an  increase 
in  maximum  lift,  rapid  variations  in  pitching  moment,  and  the  typical  lift  versus  angle 
of  attack  hysteresis  loop.  The  program  was  run  under  the  following  conditions: 

Re  =  3.45  x  io6 
Mx  =  -283 
«0  =  15° 
a,  =  10° 
k  =  .151 

For  the  first  run  the  original  Baldwin-Lomax  turbulence  model  was  used  with  an 
explicit  damping  input  of  5.  The  time  step  was  a  constant  .005  throughout,  and  data 
was  saved  beginning  at  the  mean  angle  on  the  second  oscillation  cycle.  Figure  4.4 
shows  the  lift  and  pitching  moment  coefficients  plotted  against  angle  of  attack  for  both 
theoretical  and  experimental  data.  Figure  4.5,  a  plot  of  the  pressure  drag  coefficient,  is 
included  for  completeness;  drag  plots  are  not  presented  for  the  remaining  cases,  since 
they  provided  no  additional  information  for  this  study.  Figure  4.6  contains  a  carpet 
plot  of  the  surface  pressure  coefficients  for  theoretical  data  and  a  plot  of  experimental 
data  taken  from  Reference  8.  Figures  4.7  to  4.14  are  flow-field  plots  made  at  twelve 
equal  time  intervals  during  the  cycle.  The  results  show  excellent  qualitative  agreement 
with  experimental  data  throughout  the  cycle.  The  moment  coefficient  is  consistently 
low  during  the  upstroke,  but  reasonable  quantitative  agreement  of  the  lift  coefficient  is 
obtained  until  near  18  degrees.  Just  as  in  the  steady-state  case,  the  Baldwin-Lomax 
model  gives  an  early  prediction  of  fiow  separation.  The  maximum  lift  does  not  reach 
as  high  a  value  as  experimental  results,  since  separation  occurs  at  a  lower  angle  of 
attack  in  the  theoretical  results;  as  is  the  case  for  the  experimental  data,  the  lift 
continues  to  increase  after  the  pressure  distribution,  indicated  in  the  carpet  plot,  begins 
to  break  down.  At  the  beginning  of  the  downstroke,  as  the  vortex  is  shed  off  the 
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Figure  4.4  Lift  and  Pitching  Moment  CoefFicients-Baldwin-Lomax  Model 
Mx«.283,  Re- 3.45  x  106,  a-  15#-10*cos(.3t). 
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DRAG  COEFFICIENT 


THBORZTICAL 


EXPEH  MENTAL. 


Figure  4.5  Pressure  Drag  Coefficient- Baldwin- Lomax  Model 
Mx  =  .283,  Re ™  3.45  *  106,  a»  15°-10°cos(.3t). 
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Figure  4.7  Density  Contour  Plots,  Upstroke--Baldwin-Lomax  Model 
Mx  =  .283,  Re=3.45  x  106,  a  =  15°-lO°cos(.3t) 

6.31,  9.95,  14.93,  19.99,  23.65,  25  degrees. 


Figure  4.8  Density  Contour  Plots,  Downstroke-Baldwin-Lomax  Model 
Mx-.283,  Re- 3.45  x  106,  a- 15°-10#cos(.3t) 

23.67,  20.02,  15.03,  10.03,  6.36,  5  degrees. 
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Figure  4.9  Pressure  Contour  Plots,  Upstroke-Baldvvin-Lomax  Model 
Mx=*.283,  Re  =*3.45  x  106,  a»  l5°-10°cos(.3t) 

6.31,  9.95,  14.93,  19.99,  23.65,  25  degrees. 


54 


Figure  4.10  Pressure  Contour  Plots,  Downstroke-Baldwin-Lomax  Model 
Mx-.283,  Re* 3.45  *  106,  a*  15 I0°cos(.3t) 

23.67,  20.02,  15.03,  10.03,  6.36,  5  degrees. 
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Figure  4.12  Mach  Number  Contour  Plots,  Downstroke-Baldwin- Lomax  Model 
Mx».283,  Re  =  3.45  x  106,  a  -  15*-10#cos(.3t) 

23.67,  20.02,  15.03,  10.03,  6.36,  5  degrees. 
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Figure  4.13  Stream  Function  Plots,  Upstroke-Baldwin-Lomax  Model 
Moo -.283,  Re- 3.45  x  106,  a- 15#— 10#cos(.3t) 

6.31,  9.95,  14.93,  19.99,  23.65,  25  degrees. 
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Figure  4.14  Stream  Function  Plots,  Downstroke-Baldwin-Lomax  Model 
Moo *.283,  Re=3.45  *  106,a  =  15°-10°cos(.3t) 

23.67,  20.02,  15.03,  10.03,  6.36,  5  degrees. 
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trailing  edge,  there  is  a  jump  in  all  the  experimental  coefficients.  This  is  not  duplicated 
in  the  theoretical  data.  The  flow-field  plots  reveal  the  reason  for  this  difference,  as  the 
vortex  is  not  shed  immediately  but  remains  on  the  trailing  edge  and  is  dissipated  during 
the  downstroke.  The  absence  of  strong  leading-edge  peaks  in  the  carpet  plot,  until  the 
upstroke,  shows  the  effect  of  this  discrepancy  on  the  pressure  distribution.  The 
recovery  from  separation  occurs  gradually  as  the  vortex  dissipates.  Skin  friction  values 
are  negative  over  the  upper  surface  until  the  downstroke  begins,  and  full  recovery  is 
not  observed  until  completion  of  the  downstroke. 

A  run  was  made  under  these  same  experimental  conditions  but  using  the 
modified  turbulence  model  on  the  upstroke.  It  was  intended  to  turn  off  the  modified 
model  above  19  degrees,  as  recommended  by  Sankar,  but  due  to  an  error  in  the  code,  it 
was  used  throughout  the  upstroke.  On  the  first  attempted  run,  execution  terminated  as 
the  angle  approached  19  degrees.  The  values  of  Aq*  had  begun  growing  rapidly  near 
the  leading  edge,  and  when  a  negative  value  was  computed  for  the  density,  a  math 
error  occurred  causing  the  computer  program  to  terminate  abnormally.  The  surface 
pressure  distribution  had  been  showing  unusually  high  leading-edge  peaks,  but  the 
integrated  coefficient  values  were  close  to  experimental  values  when  termination 
occurred.  Reducing  the  time  step  in  half  to  0.0025  at  15  degrees  allowed  the  solution 
to  proceed  to  above  23  degrees,  but  execution  was  then  terminated  in  the  same 
manner.  The  calculations  were  restarted  at  this  time  using  a  time  step  of  .005  and 
distance  of  the  first  point  off  the  wall  of  .00005,  but  this  time  with  the  viscosity  input, 
WW  -  10.  Reduced,  but  very  strong  leading-edge  suction  peaks  still  appeared,  the 
integrated  coefficients  were  not  as  large,  and  a  completed  solution  was  obtained 
without  reducing  the  time  step.  Data  was  again  saved  after  one  and  one-fourth 
oscillation  cycles. 

The  run  using  the  modified  turbulence  model  shows  marked  differences  from  the 
run  with  the  Baldwin-Lomax  model.  The  lift,  drag,  and  moment  coefficients  are 
plotted  along  with  experimental  values  in  Figure  4.15  Now  the  rapid  drop  in  lift  is 
delayed  too  long-until  the  maximum  angle  is  passed.  Maximum  lift  is  still 
underpredicted,  but  the  higher  viscosity  value  might  be  expected  to  affect  accuracy. 
The  moment  coefficient  is  again  lower  than  experimental  during  the  upstroke  until  the 
vortex  movement  that  is  not  predicted  by  the  code.  Note,  however,  the  agreement  with 
experiment  once  the  "theoretical"  vortex  has  been  generated.  The  carpet  plot,  Figure 
4.16,  reveals  interesting  differences  from  Figure  4.6.  Until  past  the  mean  angle  of  15 
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degrees,  the  pressure  distributions  are  nearly  identical.  In  fact,  the  Baldwin-Lomax 
model  has  a  slightly  stronger  suction  peak  (-7.7  to  -7.5  for  the  modified  model).  After 
this  point,  however,  the  Baldwin-Lomax  model  shows  the  effect  of  the  premature  flow 
separation  while  the  suction  peaks  for  the  modified  model  grow  abnormally  high.  Note 
the  different  scales  used  in  the  experimental  and  theoretical  carpet  plots.  Although  the 
integrated  lift  coefficient  begins  to  drop  slightly  nearing  the  maximum  angle,  the 
pressure  distribution  is  still  apparently  undisturbed  at  25  degrees,  with  some  weakening 
of  the  suction  peak.  The  situation  changes  abruptly  as  the  downstroke  begins  and  the 
modified  turbulence  model  is  replaced  by  the  original  Baldwin-Lomax  model.  The 
suction  peak  finally  breaks  down,  and  by  17.6  degrees  on  the  downstroke,  the 
distribution  is  again  very  similar  to  that  with  the  Baldwin-Lomax  model.  During  the 
upstroke,  the  first  negative  surface  friction  values  appear  near  the  trailing  edge  at  21.6 
degrees,  and  by  22.9  degrees,  the  negative  values  extend  over  the  upper  surface.  This  is 
the  only  clear  indication  of  the  residual  stresses  that  are  so  abruptly  relieved  when  the 
airfoil  pitches  down. 

The  flow-field  plots  are  similar  in  many  respects  to  those  for  the  Baldwin-Lomax 
model  except  for  the  timing  of  the  vortex  formation.  Figures  4.17  and  4.18  are  flow- 
field  plots  at  the  maximum  angle  and  at  23.67  degrees  on  the  downstroke.  At  the 
maximum  angle  there  is  still  no  strong  vortex  formation,  although  a  narrow 
recirculating  region  is  indicated  near  the  trailing  edge.  The  rapid  development  and 
shedding  of  the  vortex  occurs  during  the  downstroke.  By  23.62  degrees,  the  vortex  has 
progressed  well  along  the  upper  surface.  The  sequence  of  events  is  reasonable, 
although  the  timing  is  not. 

The  program  was  corrected  and  run  again  under  th*  same  experimental 
conditions  with  two  approaches  to  the  problem  of  maintaining  stability  at  the  higher 
angles.  First,  the  solution  was  marched  from  the  mean  angle  (using  the  solution  saved 
in  the  previous  run  at  15  degrees)  to  the  maximum  angle  with  the  time  step  left  at  a 
constant  .005  and  the  viscosity  increased  to  10.  Then  the  run  was  repeated  with  the 
viscosity  input  left  at  5  and  the  time  step  decreased  to  .0025  (as  had  been  attempted 
earlier  with  the  faulty  turbulence  model).  This  latter  run  was  continued  into  the 
downstroke  using  a  time  step  of  .005,  and  the  integrated  lift  and  pitching  moment 
coefficients  are  plotted  in  Figure  4.19.  The  results  are  improved  over  those  with  the 
Baldwin-Lomax  model  (Fig.  4.4).  The  maximum  lift  of  2.1  matches  experiment  and  the 
lift  curve  slope  is  sustained  longer  into  the  upstroke.  The  negative  moment  peak  is 
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ANGLE  OF  ATTACK 


Figure  4.15  Lift  and  Pitching  Moment  Coefficients- Modified  Model  until  25° 
Mx-  .233,  Re- 3.45  x  106,  a- 15°-10#cos(.3t). 
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Figure  4.17  Density  (top)  and  Pressure  (bottom)  Plots- Modified  Model  until  25* 
M so*. 283,  Re  - 3.45  *  106,  a- 15*-IO*cos(.3t) 

25°  (1),  23.67*  (r). 
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ANGLE  OF  ATTACK 


Figure  4.19  Lift  and  Pitching  Moment  Coeflicients-Modified  Model  until  19 
M 30  =  . 283,  Re-3.45  x  106,  a-  15°-I0°cos(.3t). 


sharper.  With  the  time  step  set  at  .005  and  higher  damping  (WW  -  10),  results  were 
not  as  good  guantitatively.  Peak  values  were  reached  at  near  the  same  points,  but  lift 
reached  only  1.94,  and  peak  moment  was  badly  underpredicted  at  -2.5.  For  a  final 
comparison,  the  solution  was  restarted  again  for  a  short  run  from  15  degrees,  but  this 
time  with  the  modified  turbulence  model  turned  off  immediately.  In  all  cases,  a  second 
leading-edge  suction  peak, indicative  of  vortex  form,  tion,  was  evident  in  the  first  output 
after  the  model  was  switched  off.  It  seems  from  these  results  that  the  use  of  a  smaller 
time  step  during  the  high  angle  portion  of  the  upstroke  is  worth  the  added 
computational  expense. 

C.  COMPARISON  WITH  ATTACHED  FLOW  EXPERIMENTAL  DATA 

The  program  was  next  run  under  the  experimental  conditions  of  Frame  9118  of 
Reference  8.  These  conditions  were  chosen  because,  in  contrast  to  the  previous 
experimental  results,  flow  remained  attached  throughout  the  oscillation  cycle. 
Sufficiently  lowering  the  reduced  frequency  (as  in  Frame  9110  of  Reference  8)  would 
result  in  separated  flow.  The  program  was  run  under  the  following  conditions: 

Re  *  2.45  x  io6 
Mqq  -  .184 

«o-  8# 

«!  -  10° 
k  -  .199 

Again  both  turbulence  models  were  tried,  but  for  these  conditions  it  was  possible  to 
complete  both  runs  with  a  time  step  of  .005  and  explicit  viscosity  coefficient  of  5. 
Output  was  again  obtained  after  one  and  one-fourth  cycles. 

As  was  expected,  the  Baldwin-Lomax  model  gave  a  separated  flow  prediction 
during  the  upstroke.  The  first  negative  surface  friction  values  appear  near  the  trailing 
edge  at  16.3  degrees,  and  near  17  degrees  a  vortex  forms  at  the  leading  edge.  The  lift 
and  moment  plots  of  Figure  4.20  show  this  discrepancy  from  experimental  results.  The 
lift  falls  below  experimental  values  at  high  angles,  but  again,  the  action  of  the  vortex 
sustains  the  lift  coefficient  and  even  causes  the  slope  of  the  lift  curve  to  increase 
slightly  before  the  downstroke  begins.  The  downstroke  show's  similar  behavior  to  the 
dynamic  stall  case  investigated  previously,  with  a  sharp  drop  in  lift,  which  remains  at  a 
low  value  until  the  upstroke  begins.  The  moment  coefficient  is  again  low  during  the 
upstroke  and  oscillates  as  the  vortex  propagates  along  the  airfoil  surface.  The  carpet 
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Figure  4.20  Lift  and  Pitching  Moment  Coefficients-- Baldwin- Lomax  Model 
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Moo*. 184,  Re* 2.45  x  106,  a  =  8®— 10®cos(.4t). 

plot,  Figure  4.21,  shows  the  familiar  ripple  pattern  of  vortex  movement.  The  flow  field 
at  the  top  of  the  cycle  (18  degrees)  is  illustrated  in  Figure  4.22. 

Based  on  the  previous  results  with  the  modified  turbulence  model,  it  was  expected 
to  perform  better  under  this  set  of  conditions  than  the  Baldwin- Lomax  model.  The 
flow  does  in  fact  remain  attached  during  the  upstroke  but  it  separates  as  soon  as  the 
downstroke  begins,  as  it  did  in  the  dynamic  stall  case  when  the  modified  model  was 
used  throughout  the  upstroke.  The  lift  and  drag  curves  of  Figure  4.23  are  very  similar 
in  appearance  to  those  of  Figure  4.20  except  for  the  loops  formed  at  the  maximum 
angle  by  the  rapid  vortex  formation  on  the  downstroke.  The  carpet  plot,  Figure  4.24, 
agrees  quite  well  with  experiment  until  the  maximum  angle.  The  flow  field  is  shown  by 
Figure  4.25  to  be  very  similar  at  16.64  degrees  on  the  downstroke  to  that  for  the 
Baldwin-Lomax  model  at  the  maximum  angle  (Fig.  4.22). 

D.  SIMULATIONS  UNDER  PROPOSED  EXPERIMENTAL  CONDITIONS 

A  series  of  dynamic  stall  wind  tunnel  experiments  is  presently  being  planned  at 
the  Ames  Research  Center  Fluid  Mechanics  Laboratory.  These  experiments,  w’hich 
will  include  investigation  of  compressibility  effects,  are  to  be  conducted  in  a  Reynolds 
number  range  lower  by  a  factor  of  ten  than  those  reported  in  Reference  8.  The  last 
two  simulations  attempted,  therefore,  were  made  at  a  Reynolds  number  of  345,000.  In 
order  to  have  some  basis  for  comparison,  the  mean  angle,  amplitude  of  oscillation,  and 
reduced  frequency  chosen  were  the  same  as  for  the  deep  dynamic  stall,  high  Reynolds 

number  runs.  The  Mach  numbers  used  were  .284  and  .5.  The  conditions  then  were  the 
following: 

Re  =  345,000 

Mqo  -  .284  (Run  1),  .5  (Run  2) 

«0  =  15® 

<*!  -  10® 
k  =  .151 

The  conditions  are  within  the  range  of  the  planned  experiments.  The  Baldwin-Lomax 
model  w’as  used  for  both  runs.  The  first  run  was  conducted  with  a  constant  time  step 
of  .005  and  viscosity  coefficient  input  of  5.  At  the  higher  Mach  number  of  .5,  the  time 
step  was  set  to  .0025.  An  attempt  to  start  the  run  with  a  .005  time  step  resulted  in 
abnormal  termination  after  fewer  than  30  steps  due  to  the  calculation  of  negative 
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Figure  4.22  Flow-field  Plots--Baldwin-Lomax  Model 
Mx-.l84,  Re- 2.45  x  106,  a=  8°-10°cos(.4t) 
Top-Density  (1),  Pressure  (r)  Bottom-Mach  No.  (1),  Stream  Function  (r). 
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Figure  4.23  Lift  and  Pitching  Moment  Coefficients-- Modified  Turbulence  Model 
Mx-.184,  Re- 2.45 x  106,  a- 8#-10#cos(.4t). 


72 


Figure  4.24  S 
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density  near  the  leading  edge,  just  as  had  happened  during  the  earlier  high  Reynolds 
number  run  with  the  modified  turbulence  model.  With  the  smaller  time  step,  little 
improvement  was  realized  until  the  artificial  viscosity  input  was  increased.  Fewer  than 
fifty  steps  were  completed  with  the  explicit  viscosity  at  5.  When  this  was  doubled  the 
solution  reached  nearly  8  degrees  angle  of  attack  and  completed  2170  steps  before  it 
again  "blew  up".  In  each  instance  in  which  this  abnormal  termination  took  place,  the 
output  values  all  appeared  normal  until  the  values  of  the  residuals  began  growing 
shortly  before  termination.  Increasing  the  explicit  viscosity  coefficient  to  20  permitted 
a  complete  solution  to  be  obtained. 

The  results  for  these  two  runs  were  qualitatively  similar  to  those  for  the  higher 
Reynolds  number.  The  lift  and  moment  coefficients  plotted  in  Figure  4.26  show  values 
close  to  those  of  the  earlier  results.  The  slopes  of  the  lift  curves  are  somewhat  steeper 
than  the  previous  results,  and  the  lift  continues  to  drop  (and  the  moment  coefficient  to 
rise)  farther  into  the  downstroke.  Figure  4.27  compares  the  pressure  distribution  at  the 
lower  Mach  number  with  earlier  results  for  the  high  Mach  number.  The  suction  peaks 
at  the  low  Reynolds  number  are  not  as  strong  and  begin  breaking  down  slightly  earlier. 
The  flow-field  plots  showed  a  similar  series  of  events  to  those  seen  previously  and  are 
not  presented  here. 

The  results  obtained  at  the  higher  Mach  number  have  two  notable  features: 
smoother  contours  and  lower  peak  values.  At  least  some  of  this  smoothing  must  be 
attributed  to  the  much  higher  viscosity  input  used  for  this  case  (20  versus  5).  The  lift 
and  moment  plots  are  given  in  Figure  4.28  and  the  carpet  plot  for  this  case  along  with 
that  at  the  lower  Mach  number  are  shown  for  comparison  in  Figure  4.29.  The  shape 
of  Figure  4.28  curves  is  more  rounded  than  was  Figure  4.26,  and  maximum  lift  is 
reduced.  The  pressure  coefficient  suction  peaks  (Fig.  4.29)  are  also  smoother,  lower, 
and  begin  breaking  down  sooner  at  the  higher  Mach  number.  To  further  investigate 
the  effect  of  artificial  viscosity,  a  partial  run  was  made  with  the  explicit  damping 
coefficient  input  as  55.  The  maximum  lift  coefficient  obtained  was  1.3  at  17.6  degrees 
angle  of  attack  compared  to  1.6  at  18.2  degrees  with  the  viscosity  at  20.  The  maximum 
value  was  attained  at  nearly  the  same  angle  (output  was  not  taken  at  exactly  the  same 
points  for  this  second  run).  However,  the  first  irregularities  in  the  leading-edge 
pressure  coefficient  were  delayed  by  this  increase  in  viscosity  from  about  12.3  degrees 
angle  of  attack  to  near  15  degrees.  Thus,  the  computer  results  do  seem  to  be 
predicting  an  earlier  vortex  formation  (and  breakdown  in  the  surface  pressure 
coefficient)  at  the  higher  Mach  number. 
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Figure  4.26  Lift  and  Pitching  Moment  CoefTIcients-Baldwin-Lomax  Model 
M x -.284,  Re-. 345  x  I06,  a- 15°— 10#co$(.3t). 
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Figure  4.28  Lift  and  Pitching  Moment  CoefFicients-Baldwin- Lomax  Model 
M  oo  -  .5,  Re  -  .345  x  lo6,  a  -  15*-10*cos(.3t). 
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V.  CONCLUDING  REMARKS 


The  computer  program  employed  for  this  study  has  shown  the  capability  to 
model  the  basic  events  of  the  dynamic  stall  process.  The  quality  of  the  results  is,  not 
surprisingly,  highest  at  moderate  angles--two  to  three  degrees  above  static  stall  angle  in 
examples  studied.  Results  at  higher  angles  are  strongly  dependent  on  turbulent 
viscosity  model.  The  downstroke  is  not  as  well  approximated,  although  there  is 
qualitative  agreement  of  the  integrated  coefficients.  The  program  appears  quite  robust 
at  Mach  numbers  of  .3  and  below  using  the  Baldwin-Lomax  turbulence  model.  With 
the  modification  to  the  turbulence  model  introduced  in  this  implementation,  better 
results  are  obtainable  at  the  higher  angles,  but  the  program  is  less  stable.  Adequate 
convergence  was  generally  achieved  quickly  enough  that  data  could  be  saved  beginning 
at  the  mean  angle  on  the  first  oscillation  cycle. 

The  examples  considered  show  some  of  the  limitations  of  current  Navier-Stokes 
solvers.  By  manipulating  artificial  viscosity  and  turbulence  modelling,  it  is  possible  to 
manage  calculations  to  match,  in  some  degree,  desired  results.  For  example,  the 
attached  flow  case  considered  might  be  brought  into  closer  agreement  with  experiment 
by  using  the  modified  turbulence  model  after  the  beginning  of  the  downstroke.  For  the 
dynamic  stall  case,  however,  this  might  lead  to  premature  flow  reattachment  and  less 
accurate  results.  The  simple  modification  introduced,  used  as  recommended  during  the 
upstroke,  shows  encouraging  results.  At  no  significant  expense,  modelling  near  and  a 
few  degrees  above  static  stall  angles  appears  improved  without  affecting  results  at 
lower  angles.  Further  testing  needs  to  be  done  to  determine  the  effectiveness  of  the 
change  under  different  dynamic  conditions,  and  the  use  of  other  models  is  certainly 
worth  consideration.  King  and  Johnson  have  reported  favorable  results  with  a 
nonequilibrium  turbulence  model  using  an  ordinary  differential  equation  to  model 
streamwise  shear  stress  development,  and  an  eddy  viscosity  distribution  across  the 
boundary  layer  based  on  the  maximum  value  from  this  equation.  Compared  to  the 
Baldwin-Lomax  and  Cebeci-Smith  models,  it  showed  little  difference  in  attached  flow 
calculations  with  a  Navier-Stokes  solver  using  a  Beam-Warming  scheme.  When 
applied  to  separated  flows  and  flows  w’ith  strong  viscous  effects,  it  gave  superior 
results,  though  at  some  computational  expense  in  its  present  form.  [Refs.  18,20] 
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The  use  of  other  grids  or  variations  on  the  present  one  would  be  worth 
considering  for  the  purpose  of  developing  comparisons  for  the  proposed  experiments  at 
the  Ames  Research  Center  Fluid  Mechanics  Laboratory.  Obtaining  a  solution  at  the  .5 
Mach  number  was  difficult  and  required  a  large  increase  in  the  artificial  viscosity. 
Some  of  the  problems  might  be  alleviated  by  closer  spacing  in  the  leading-edge  area 
where  the  instability  arose.  The  present  grid  generation  scheme  might  be  tried  with 
more  points  selected  in  the  normal  direction.  Events  at  the  trailing  edge  are  not  now 
modelled  with  complete  accuracy.  The  vortex  moves  to  the  trailing  edge  and  diffuses 
there.  This  behavior  may  be  influenced  by  the  level  of  artificial  viscosity  (which  must 
be  increased  to  compensate  for  grid  coarseness),  the  downstream  boundary,  location  of 
the  grid  cut,  and  relative  coarseness  of  the  grid  at  the  trailing  edge.  The  present  grid 
provides  adequate  results  for  a  large  range  of  applications  with  a  reasonable  number  of 
grid  points,  however. 

The  solver  now  assumes  a  fully  turbulent  boundary  layer  (or  fully  laminar  if  the 
Reynolds  number  is  set  to  zero).  The  proposed  experiments  may  include  significant 
transition  point  effects.  The  transition  could  be  simulated  by  retaining  the  laminar 
viscosity  coefficient  forward  of  a  specified  chord  location.  This  might  be  especially 
useful  for  comparison  with  tripped  boundary  layer  data.  Some  runs  were  made  with 
fully  laminar  viscosity.  The  pressure  gradients  at  the  leading  edge  could  not  be 
sustained,  and  multiple  vortices  were  formed  but  the  solution  remained  stable. 

The  cases  considered  in  this  study  represent  useful  test  conditions  for  studying 
the  effect  of  changing  various  program  segments.  The  deep  stall  case  has  the  model 
dynamic  stall  features  and  has  been  used  by  Sankar.  The  attached  flow  case  is  an 
interesting  and  a  difficult  contrasting  case.  The  forthcoming  experiments  will  provide 
the  opportunity  to  investigate  details  of  the  dynamic  stall  prdcess,  including 
compressibility  effects,  in  much  greater  depth.  If  the  Navier-Stokes  solver  is  not  yet  a 
completely  adequate  predictive  tool,  the  effort  to  make  it  so  can,  of  itself,  provide 
insight  into  the  nature  of  unsteady  flows. 
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APPENDIX  A 

SANKAR  NAVIER-STOKES  SOLVER 


JO0.JN-.T-. 

ACCOUNT.  AO,  US-,  UPW-, 

• . 

CFT.ON-A.OfF-S. 

• !  ASSIGN.  ON-NEWS  LN.A-FTM. 

•  .  ASSIGN. 0N-NEWCP.A-FT18. 

• .  ASSIGN .  ON-XYZ .  A-FT1 1 . 

*. ASSIGN, 0N-Q.A-FT12. 

•  .IS  THIS  RUN  TO  START  FROM  A  STORED  SOLUTION? 

•  .  ACCESS .  ON-OLDSLN ,  PON- ,  I D- . 

•  .00  YOU  WANT  TO  ACCESS  EXISTING  PRESSURE  COEFFICIENT  OATA? 

•  .  ACCESS .  ON-OLOCP .  PON- ,  ID- . 

•  !  ASSIGN. DN-OLOSLN.A-FTB7. 

•  .ASSIGN.  ON-OIDCP  .A-FT17. 

LOR .  MAP-PART .  SET-ZERO . 

•  DO  YOU  RANT  TO  SAVE  THE  CURRENT  SOLUTION? 

• .  SAVE .  DM-NENSLN .  PON- .  ID- . 

•  .DO  YOU  WANT  TO  SAVE  PRESSURE  COEFFICIENT  DATA? 
•.(PREVIOUS  FILES  SHOULD  BE  PURGED  FREQUENTLY). 

•  .  SAVE .  DN-NEWCP .  POIN ,  ID- . 


•  DO  YOU  WANT  TO  CREATE  PL0T30  FILES? 

•  .ACCESS.  ON-SENDVAX,  POfNSENDVAX.  ID-STTROM. 

•  .  SENDVAX ,  DfNXYZ ,  VDfN . 

•  .  SENDVAX.  DN-Q,  VON- . 

•  DO  YOU  WANT  PL0T3D  FILES  FOR  TROUBLESHOOTING  IF  PROGRAM  CRASHES? 
•EXIT. 

•  .ACCESS. OINSENDVAX.PON-SENDVAX,  ID-STTROM. 

•  .  SENDVAX.  DN-XYZ.  VON- . 

•. SENDVAX. DN-O.VON-. 

•  !00  YOU  WANT  TO  USE  JOB  CHAINING? 

•  .  FETCH,  WN,  TEXT-. 

•.REWIND, ON-. 

•.SUBMIT. ON-. 

/EOF 

PROGRAM  MAIN 

•  PROGRAM  MAIN(  INPUT  .OUTPUT ,  TAPES- INPUT .  TAPE6-0UTPUT) 
COMJON/SURF/PSUR(  1  •  1 ) 

COMJON/F I X/OMEGA 
COAMON/MUTUR/CMU  (161 .41) 

CO»MON/SKINCF/CF(161) 

comjon/gridiA(i«i.ai).2(i«i.4i) 

COMJON/PAR/GAN4A.REYREF.  ALFA,  ALFA  1 ,  REDFRE ,  AMINF ,  ALFAI 
COMJON/OC*  IO/DT .  I  MAX ,  KMAX ,  I  TEL .  I TEU 
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C0MM0N/GRID/YAC08(  161 .41 ) 

COUCN/DAkP/WW ,  WW I 

COIMON/FLOW/Q1 (161 .41),  02(161, 41), 03(161, 41)  ,04(161 ,41) 

•' ALPHA<®«)  .CPTH(97.*«)  ,XTH(97) 
COMON/LOG I C/RSTRT ,  P I TCH ,  PLUNGE 
LOCICAL  RSTRT, PI TCH, PLUNGE 
C 

Cl  1 1 1 1  INDICATES  COMMENTS  OR  CHANGES  ADDED  BY  J.F.  VALDES.  DEC  1986 
C*»*  PROGRAM  SOLVES  TWD-OIMENSIONAL  FLOUT  PAST  ARBITRARY 
GEOMETRIES  USING  AO  I  PROCEDURE. 


C*** 

c 


c 

c 

c 

c 

c 

Cl  I 
Cl  I 
Cl  I 
CM 
Cl  I 
Cl  I 

c 

c**« 

c 


TAPES  -  FILE  CONTAINING  INPUT  OATA 
TAPES  -  OUTPUT 

T^pE8-FILE  THAT  SAVES  THE  FLOW  FIELD  AT  THE  END  OF  A  RUN 
IF  THE  CURRENT  RUN  IS  A  RESTART  OF  A  PREVIOUS  RUN,  THEN 
T0  REA0  THE  FL0"  FIEL0  into  MEMORY 
T  ?i8JirU!P.I2,i$ClMJLATE  PRESSURE  COEFFICIENT  OATA  FOR  A  CYCLE 

"n®1  C0MPLETE  by  PROGRAM  PLOTNSE.SRC 
fJLJf  USED  JO  READ  EXISTING  PRESSURE  COEFFICIENT  DATA  FOR 
1 1 1  UPDATE  OUR  INC  THE  CURRENT  PROGRAM  RUN 
I  I ITAPE11  IS  USED  FOR  PL0T30  XYZ-FILE  OUTPUT 
I  I  ITAPE12  IS  USED  FOR  PL0T3D  O-FILE  OUTPUT 

READ  INPUT  OATA 

PI  -  ATAN( 1 . )*4. 

READ  jS.lj  TITLE 


READ 


-  IMAX.KMAX.OT.WW.ALFA.ALFAI  .ALFAI  .REDFRE.AMINF 

Cl  I  I  I  I  SET  ALFAD  FOR  STEADY  STATE  PL0T3D  OUTPUT 
ALFAD  -  ALFA 


C 

C 

C 

C 


C 

C 

C 


FORCE  DT  TO  BE  EQUAL  TO  UNITY  FOR  STEADY  FLOW  PROBLEMS 
THIS  INVOKES  THE  RELAXATION  MOOE 

IF(REDFRE. LE.6.601 )  DT  -  1,6 
READ  (S.2661) 

s,eps  io  eE «« *  ™is «« 

NSTP  -  FNSTP 
READ  (3.2661) 

REYREF-  REYNOLDS  NUMBER  IN  MILLIONS 

^!N  -  DISTANCE  OF  FIRST  POINT  OFF  THE  WALL 

FOR  REYNOLDS  NUflERS  UPTO  3  MILLION  USE  6.66665 

READ  (S.2221 )  REYREF. DNMIN 

REYREF  -  REYREF  •  1 . E+66 


(5.2661) 
(3.2221)  TSTART 


JSTART  -  TIME  THAT  THE  CALCULATIONS  HAVE  BEEN  ADVANCED 

OBTAINED  FrS’taPeT'  ^  TSTA*T  IS  NECATIVE  ™,S  VALUE  IS 
READ 

READ  _ 

C 1 1 1 1 1 FULOUT  IS  A  FLAG  FOR  GENERATING  PLOTTING  FILES.  WHEN  NEGATIVE 
CIIIIINO  DATA  IS  GENERATED.  ZERO  IS  USED  TO  START  AND  A  POSITIVE 

C' '' ' '»r«W«T?«2??TINUE  GENERAT,NC  fU>-L  OUTPUT.  ADOED  FEATURE( JFV)  . 
(9 1 2M i  i 

READ  (3,2221)  FULOUT 
2221  F0RMAT(2F16.4) 

READ  (5,2661) 

2661  FORMAT(IX) 

READ  (5.2666)  RSTRT, PI TCH, PLUNGE 
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ooo 


2999  F0RMAT(4L5) 

C  NEGATIVE  REYREF  MEANS  INVISCIO  FLOW 
C 

C*««  PRINT  OUT  THE  INPUT  DATA 
C 

WRITE  (9,4)  TITLE 

WRITE  (6,3)  IMAX.KMAX.DT.WW.ALFA.ALFAI  .ALFAI  .REDFRE.AM1NF 
IF(REYREF.GT.9. )  WRITE  (6.3799)  REYREF 
GAI64A-1.4 
C 

Cl  1 1 1 (GENERATE  COMPUTATIONAL  GRIO 
C 

CALL  AIRFOL(  IMAX.KMAX,  1TEL,  ITEU) 

IF(REYREF.GT.9.)CALL  CLUSTR(DNMIN) 

WRITE  (6,3992) 

C 

C«««  STARTING  CONDITIONS. 

C***  DENSITY  NORMALISED  WITH  RESPECT  TO  ROINF 

C**«  VELOCITIES  NORMALISED  WITH  RESPECT  TO  AINF 

(*•••  TOTAL  ENERCY  NORMALISED  WITH  RESPECT  TO  (ROINF*AINF*AINF) 

C 

TOTE^AMINF«AMINF«9.54>1  ./(GAM4A*(GAM4A-1 . )) 

ALFA  -  ALFA  .  ATAN(1 . )  /  43. 

ALMEAN  -  ALFA 

ALFAI  -  ALFAI  •  ATAN( 1 . )  /  43. 

ALFAI  •  ALFAI  •  ATAN( 1 . )  /  45. 

ALFAS  -  ALMEAN  -  ALFAI 
UINF  -  AMINF  •  COS(ALFA) 

VINF  -  AMINF  •  SIN(ALFA) 


DO  7  1-1 . 1  MAX 

DO  7 

(-1.KMAX 

01 

I.K, 

-1 . 

02 

I.K 

-UINF 

03 

I.K 

-VINF 

04  ( 

I.K) 

-TOTEN 

7  CONTINUE 

IF(RSTRT)  READ  (7)  TIME, 01 ,02,03.04 
IF(TSTART.GE.9. )  TIME  -  TSTART 
IF( .NOT. (RSTRT))  TIME  -  6. 


1 1 1 1  (READ  STORED  PRESSURE  COEFFICIENTS.  STATEMENTS  ADO  ED  BY  JFV. 


IF  (FULOUT  .GT.  9.)  THEN 
READ  (17,1)  TY  LE 

READ  (17)  AMPLTD.STMEAN, CMS. XTH, ALPHA. CPTH 
END  IF 

Clll  I  (BEGIN  PRESSURE  COEFFICIENT  DATA  FILE.  STATEMENTS  ADO  ED  BY  JFV 
IF  (FULOUT  .EQ.  9.)  THEN 

AMPLTD  •  ALFAI  •  43.  /  ATAN(1 . ) 

STMEAN  •  ALFA  •  43.  /  ATAN(  1 . ) 

OMS  -  REDFRE 
TEDGE  -  X(ITEL+48,1) 

00  IS  1-1,97 

XTH(I)  -  X(I+ITEL-1.1)  -  TEDGE 
IS  CONTINUE 
END  IF 


Cl  1 1 1  (DETERMINE  NUMBER  OF  TIME  STEPS  OF  CURRENT  SIZE  IN  ONE  CYCLE 
Clll  1 1  SCHEME  FOR  PLOTTING  DATA  A  NO  PROGRAM  EXIT  ADOED  BY  JFV 
IF  (PITCH)  CYCLE  -  PI  /  (REOFRE  •  AMINF  •  DT) 

C  •••• 
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oooo  o  o  o 


c 

CALL  METRIC 
00  It##  ITN-1.NSTP 
C 

Cl  III  (CREATE  PLOT 30  FILES  FOR  INVESTIGATING  INSTABILITY  ADOED  BY  JFV 
Cl  1 1 1 1  THIS  BLOCK  UPOATES  PL0T30  OUTPUT  FILE  AFTER 
Clllll  EACH  STEP.  USE  IN  CONJUNCTION  WITH  EXIT  JCL. 

Cl  II  1 1  ••••  NOTE  -  USE  ONLY  FOR  TROUBLESHOOTING  •••• 

C 

C  REWIND(II) 

C  WRITE(II)  IMAX.  KMAX 

C  WRITE(ll)  ((X(I,J),  1-1,  IMAX),  J-1.KMAX). 

C  1  ((Z(I.J),  1-1, IMAX),  J-1.KMAX) 

C  REWIN0(12) 

C  WRITE? 12)  IMAX.  KMAX 


C  WRITE(12)  AMINF,  ALFAO.  REYREF.  TIME 

C  WRITE(12) 

[Oil 

[I.J 

),  i-i.imax; 

I.  J— 1 , KMAX) 

C  1 

Q2 

I.J 

1,  1-1, IMAX 

.  J— 1 , KMAX) 

C  2 

i 

031 

1.  I-I.IMAX 

,  J-1.KMAX) 

C  3 

1 

041 

I.J 

.  I-I.IMAX] 

.  J-1.KMAX) 

TIME  -  TIME  +  OT 

Cl  1 1 1 1  ROTATE  GRID  TO  NEW  ANGLE  OR  SET  TO  CORRECT  ANGLE  FOR  RESTARTS 
IF  (PITCH)  THEN 

OMEGA  -  2.  •  REDFRE  •AMINF*SIN(REDFRE  •  2.  •  TIME  *  AMINF) 
1  *ALFA1 

ALOLO  -  ALMEAN  -  ALFA1  «  C0S(2.  •  REDFRE  •  AMINF  • 

1  (TIME  -  OT)) 

ALFA  -  ALMEAN  -  ALFA1  •  COS  (REDFRE  •  2.  •  TIME  •  AMINF) 
ALFAD  -  ALFA  •  45.  /  ATAN(1 . ) 

OALFA  -  ALFA  -  ALOLD 

IF(RSTRT .AND.  ITN.EQ.  1 )  OALFA  -  ALFA  -  (  ALMEAN  -  ALFA1  ) 
CALL  ROTGRIO(X.2.IMAXiKMAX.DALFA) 

CALL  METRIC 
END  IF 

IF  (PLUNGE)  THEN 

OMEGA  -  2.  •  REDFRE  •  AMINF 
ALFA  -  ALMEAN 
END  IF 

C 

Cl  I  I  I  I  COMPUTE  THE  SOLUTION  BY  ADI  TECHNIQUE 
CALL  SLPS( I TN, OMEGA, OALFA) 

APPLICATION  OF  BOUNDARY  CONDITIONS 

CALL  WALLBC 


PRINT  OUT  PRESSURE  AT  THE  SURFACE 

1 1 1 1  IFOR  STEADY  STATE  OUTPUT  USE  THE  FOLLOWING 

IF((ITN/1##*1BB.EQ.ITN)  .ANO.  ( .NOT. PITCH))  THEN 
WRITE  (#.19) 

WRITE  (#,33)  ITN 

CALL  CPPLOT(  ITEL.  ITEU, AMINF, X(  1 . 1 )  .2(1 . 1 )  ,PSUR) 
WRITE  (6,12)  (I. CF(I). 1-1, IMAX) 

CALL  LOAD(PSUR,CL.CO,CM,ALFAS) 

WRITE  (6,3###)  CL.  CO.  CM 
WRITE  (6,3##2) 

END  IF 

CHIIIOR.  FOR  DYNAMIC  SOLUTION  OUTPUT  AT  EQUAL  PHASE  INTERVALS 
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CIII 1 1  MOOIFIED  OUTPUT  SCHEME  ANO  PROGRAM  EXIT  SCHEME  BY  JFV. 

IF  (PITCH)  THEN 

TINCR  -  TIMC/DT  -  CYCLE/4. 

IF  (TINCR  .IT.  •.)  TINCR  -  TINCR  ♦  CYCLE 
CIII  1 1  DETERMINE  NU0ER  OF  STEPS  BETWEEN  OUTPUT 
NCPOUT  -  HINT (CYCLE/96.) 

Cl  III  IMULTIPLY  NCPOUT  BY  THE  NUMBER  OF  OUTPUT  CYCLES  DESIRED  BETWEEN 
CIII  I  (PROGRAM  EXITS.  THIS  DETERMINES  THE  INTERVAL  FOR  PLOT3D  FILES. 
NEXIT  -  24  •  NCPOUT 
ACYC  -  AMOO(TINCR, CYCLE) 

DO  IB  J-1,4 


C 

C 


IF  (A8S(T!NCR  -  CYCLED)  .LT.  B.9  .AND.  (ITN  .OT. 
NCPOUT))  ACYC  -  96.  /  (NEXIT  /  NCPOUT)  •  NEXIT 
IUE 

N8CYC  -  N1 NT (ACYC) 

IF  (NBCYC/NCPOUT  •  NCPOUT  .EO.  NBCYC)  THEN 
INDEX  -  NBCYC  /  NCPOUT 
ALPHA( INDEX)  -  ALFAD 


•  ••• 


WRITE  (6.19) 

WRITE  (6.33)  ITN 
MUTE  (6,3566)  ALFAD 

CALL  CPPLOT(  ITEL,  ITEU.AMINF  '*(1 , 1)  .2(1 . 1 )  .PSUR) 

C 

C 1 1 1 1 1  STORE  CURRENT  PRESSURE  COEFFICIENTS.  ADDED  STATEMENTS  BY  JFV. 
DO  26  J-1 .97 

CPTH(J. INDEX)  -  PSUR( J+ITEL-1 ) 

26  CONTINUE 
C 

Clllll 


C 

C 

C 


IF  (FULOUT  .CE.  6.)  KMITE(6.12)  (I ,CF(I) . 1-1 . IMAX) 
CALL  LOAD(PSUR.CL.CD.CM.ALFAS) 

WRITE  (6.3666)  CL  .  CO  .  CM 


1 1 1 1  I  FOR  AUTOMATIC  PROGRAM  EXIT  DURING  FINAL  OUTPUT  (JFV) 


IF  ((NBCYC/NEXIT  .  NEXIT  .EG.  NBCYC)  .AND.  (ITN  .GT. 
1  2  •  NCPOUT))  GO  TO  1661 

WRITE  (6.3662) 

C  •••• 

C 


END  IF 
END  IF 
C 
C 

1666  CONTINUE 
Clllll 

1661  CONTINUE 

WRITE  (8)  TIME.Q1 .02.03.04 
CIIIIISAVE  PRESSURE  COEFFICIENT  DATA  (JFV). 

IF  (FULOUT  CE.  6.)  THEN 
WRITE  (18.1)  TITLE 

MUTE  (18)  AMPLTD . STMEAN , OMS , XTH , ALPHA , CPTH 
C 

Cl  III ICREATE  PL0T30  FILES.  FEATURE  ADOED  BY  JFV. 

Cl  1 1 1  IDO  NOT  USE  MIEN  TROUBLESHOOTING  (ABOVE) 

WRITE(II)  IMAX,  KMAX 

WRITE(II)  ((X(I.J),  1-1, IMAX),  J-1. KMAX). 
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OUU  UU 


((Z(I.J),  1-1 . IMAX) .  J-1.KM4X) 
WWITCf 12)  I MAX ,  KMAX 
WRITE(12)  4MINF.  ALFAO,  REYREF.  TIME 


1 

2 

3 


WRITE 


I-1.IM4X) 
1*1 , 1  MAX) 


J-1 , KMAX 
J-1.KMAX 


1*1 , 1  MAX) ,  J-1.  KMAX 
1-1 , 1  MAX)  ,  J-1  .KMAX 


END  IF 

PRINT  OUT  VELOCITY  PROFILE 

00  4»H  1  -  1  .  IMAX  .  2 

S  -  0. 

DO  4«M  K  -  2  .  10 

*  -  s «-.5?rI^x^,  k)‘x(i  k-1»**2+(z(,  k)“2(i-k-i))**2) 

^  ■  CMUi  I  (K)  /  QT 

UTOT  -  SORT (02 ( I , K)«»2  ♦  Q3(  I  ,K)**2)/Q1 (I ,K) 

I II 1 1  IF  PRINTED  OUTPUT  IS  DESIRED  ' 

400#  CONTINUE*  (#‘2##2)  1  *  KlQ2(I*K)  •  M(I.K)  .  ED  .  UTOT , S 
STOP 

1  FORMAT (15A4) 

2  FORMAT(2I9,7F10.0) 

3w°SfAli/*221i5^I5fA5“v,2*-//'5X-5HK,-AX"'  OT-.F20.8. 

V/ '  55' ■  F£LV  *,5! ‘ 9HALFAp' ■ F2B 8 > f ** .  W4LFA1-, „ F20 .1 0 . 

4  5X * 7HREOn,E" . P20 . 8 ./. 5X . 8HAMINF-, F20.8) 

12  FORMAT 
19  FORMAT 
22  FORMAT 
33  FORMAT 
2002  FORMAT 
3000  FORMAT 
3002  FORMAT 


1H1.SX.19A4) 

8jI4.F10.4)) 

2F10.0.I9) 

9X.9MISTP-, 19) 

2I9.9E14.0) 

3X.3HCL-.F1O.4.9X.3HCO-.F10.4,3X.3HC»fKF10.4) 
VA^’W^'.IIX.’OWIAX'.IIX.-DVMAX'.IIX.-DEMAXMRX. 

3500  FORMAT 1/ . 5X , 9HALFA- , FI 0 . 4) 

3700  FORMAT (3X,7HREYREF-,F1 3. 4) 

END 

SUBROUTINE  AMAT1  (K.  IMX1  .XIX.XIZ.XIT) 

>  -?2( •«’ )  .03(1*1.41 )  .04(101.41 ) 

1!  .*1). 002(101, 41).DQ3(101.41).DQ4(101, 41) 

WCY"CP -ALFA .  ALFA1 .  KEOPWE .  AMINF .  ALFAI 
DIMENSION  XIT(101.41),XIX(181,41),XIZ(101,41) 

REAL  K0.K1.K2 
C 

C...  AMAT1  COI0*UTES  THE  COEFFICIENT  MATRIX  DE/DO  DURING  XI  SWEEP 

GM1  -  GAMJA  -  1 . 

DO  1000  I  -  2  ,  IMX1 


K0 

K1 

K2 

U 

W 

EBYR 
PH  1 2 
THETA 
4(1.1. I) 
4(1.2. I) 


-  XIT(I,K 
XIXfl.K 


-  XIZ( I ,K, 

-  02(1 ,K)  /  Q1 

-  03(1. K)  /  01 

-  Q4(I.K)  /  Q1 

-  0.9  •  GM1  • 


-  K1  •  U  ♦  K2  •  W 

-  K0 

-  K1 


r.x) 

E.K) 

,  I .«) 

U  *  U  *  *  ■  ») 
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AM. 3. 1)  -  K2 
Al  1,4.1  -0 

Al  2,1,1  I  -  K1  •  PHI2  -  U  •  THETA 

Al  2,2,1  I  -  X*  ♦  THETA  -  K1  •  (OKI  -  1.)  •  U 

Al  2,3, 1)  -  X2  •  U  -  GM1  •  K1  •  W 

A(2,4, 1  -  K1  •  OM1 

A(3.1 . t 1  -  K2  •  PHI2  -  «  •  THETA 

A(3,2, I  ■  K1  •  W  -  K2  •  GM1  *  U 

A(3,3, 1  -  X*  ♦  THETA  -  K2  •  (OM1  -  1 . )  •  W 

A(3,4, I)  -  K2  •  GM1 

Al  4, 1 , 1 1  -  THETA  *  (2.  •  PH12  -  GAMMA  •  EBYR) 

Al  4,2, I )  -  K1  •  (GAMMA  •  EBYR  -  PHI2)  -  GM1  •  U  •  THETA 

Al  4,3, 1  -  K2  •  (GAMMA  •  EBYR  -  PH12)  -  GM1  •  W  «  THETA 

A(4,4, I  -  K0  ♦  GAMMA  •  THETA 
100*  CONTINUE 
RETURN 
END 

SUBROUTINE  AMAT2(I  ,KMX1  .ZETAX.ZETAZ.ZETAT) 

COMMON/ELOR/Q1  ( 101 ,41 )  ,Q2(161 ,41 )  ,03(101 ,41 )  ,04(101 ,41 ) 
C0MM0N/PERTR/DQ1  (161 ,41 )  ,002(101 ,41)  ,003(161 ,41)  ,004(101 ,41) 
C0M40N/AM/A(4,4, 161 ) 

COMMON/PAR/CAMMA ,  REYREF .  ALFA ,  ALFA1 ,  REOFRE ,  AMI  NF .  ALFA  I 
DIMENSION  ZETAX(161 ,41)  ,ZETAZ(161 ,41)  ,2ETAT(101 ,41) 

REAL  K0.K1 ,K2 
C 

C*«*  AMAT2  COMPUTES  THE  COEFFICIENT  MATRIX  OF/DO  DURING  ETA  SWEEP 

C 


AM.1.K) 
AM.2.K) 
AI1.3.K  I 
A1.4.K 
A(2.1.K  i 
AI  2.2.K  i 
AI2.3.K  I 
A(2.4.K  i 
Al  3.1.K 
AI3.2.K  i 
AI  3.3.K  i 
Ai 3.4.K 
Al  4,1,K  i 
Ai  4.2.K  i 
A(4.3.K) 
A(4.4.K) 
10*0  CONTINUE 


K* 

K1 

K2 

* 

XI  •  PHI2  -U  •  THETA 
X*  ♦  THETA  -  XI  •  (GM1-1 . )  •  U 
K2  •  U  -  GM1  •  XI  •  W 
XI  •  0M1 

X2  •  PH I 2  -  W  •  THETA 
XI  •  W  -  K2  •  GM1  •  U 
X*  ♦  THETA  -X2  •  (0M1-1.)  •  W 
K2  •  GM1 

THETA  •  (2.  •  PH  1 2  -  GAMMA  •  EBYR) 

XI  •  (GAMMA  •  EBYR  -  PHI2)  -  GM1  •  U  •  THETA 

X2  •  (GAMMA  •  EBYR  -  PHI2)  -  GM1  •  W  •  THETA 

X*  4  GAMMA  •  THETA 


RETURN 


ENO 


SUBROUTINE  SLPS(ITN, OMEGA, OALFA) 

COMMON/FLOW/01  ( 101 , 41 )  ,02(  101 , 41 )  ,03(  101 .41 )  ,04(  161 .41 ) 
COMMON/PERTR/DOI  ( 1 01 . 41 ) .  D02(  1 01 , 4 1 ) . D03(  1 01 . 41 ) , 004(  1 01 . 41 ) 
COMMON/ AM/ A  (4, 4, 101 ) 
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oo  OOO  O  O  OO 


1^CC4N1«7I2'i?0^^’** ,ei  *41)  1«1 .41)  ,Ee(^.^.  1«1 .41) 

* ALFA  * ALrA  1  • REDFRE  • AMI NF  •  A*-FA  I 
CMOH/OG*  IO/OT ,  IMAX ,  KMAX ,  I  TEL ,  I TEU 
CCM40N/CRIDAAC08(1«1,41) 

COM40N/DAMP /WW ,  WWI 

REAL  IM  ‘ 

DIMENSION  0ELTAT(161 ,41) 

0 

£•••  SUBROUTINE  SIPS  DOES  THE  BULK  OF  THE  PORK  FOR  THE  ADI  ai  nmirui 
C*«*  IT  CALLS  FLUX  AND  COMPUTES  RIGHT  HAND  SIDE  DURING  THE  TWO  SWEEPS 
C...  ASSEMBLES  THE  COEFFICIENT  MARICES^S  II^K  T  ANO  EXP?i?ff  ’ 

§:::  SOLUTION*0"  ^  CAlLS  TH€  TRIOIAOONAL  SOLVER  TO  OBTAIN  THE  FINAL 

C,H,,2EJ  VALUE  OF  IMPLICIT  OAMPING  COEFFICIENT 
RRI  ■  2B.B  •  WW 
IM1  -  IMAX  -  1 
I M2  -  IMAX  -  2 
KM1  •  KMAX  -  1 
KM2  -  KMAX  -  2 
IF(ITN.EO.I)  THEN 

. . steaov  sT4Tt  “icuu7ioHs 

DO  777  K  -  2  .  KMAX  -  1 
DO  777  I  -  2  .  IMAX  -  1 

777  CONTINuJ’K)  “  *  5  f  (1-  +  «RT(ABS(YAC08(I.K)))) 

ELSE 

DO  778  K  -  2  .  KMAX  -  1 
DO  77B  I  -  2  .  IMAX  -  1 

778  DELTAT(I.K)  -  1.0 
END  IF 

ENO  IF 
C 

C::.‘  S3! . SoISiSSTi§2  TERMS  ARE  C0NSTRUCTED  AND  STORED  IN  THE  ARRAYS  DQ1 . 


•  •• 


CALL  OISSIP 

CALLRRKI  (OMEGA) IDE  AT  KN0WN  T!ME  LEVEL  IS  M0*  COATED  AND  ADDED 

IF  VISCOUS  FLOW  IS  COMPUTED  THE  VISCOUS  TERMS  ARE  ADDED  TO  DQ1  ETC.  HERE. 
fFX2^2EF  GT  a* J  CALL  STRESS(ITN.DALFA) 


DTH  -  DT  •  0.S 

DTW  -  DT  •  WWI 

DO  3  K  •  2  ,  KM1 

CALL  AMAT1(K, IMAX-1 ,XIX,XIZ,XIT] 

00  4  II  -1,4 

DO  4  12  -1.4 

DO  S  I  -  2  ,  IMAX  -  1 

KCIJ.I2.I-1.K1  -  AfI1.I2.lii; 

00(11.12, 1-1. K)  -  -  A( 1 1 , 12, 1-i l 

CONTINUE  1 

CONTINUE 


•  DTH  •  DELTAT? I ,K) 

•  DTH  •  DELTAT ( I ,K) 
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c 

C*««  IMPLICIT  DAMPING  AOOCD  HERE. 
C 


C 

C 

C 


C 

C 

C 


00  6  II 
DO  •  I 
DTI 

00(11,11. 1-1 .K 
EEC  11.11. 1-1 ,K 
MM(I1 , II , 1-1 ,K 
CONTINUE 
00  DM  I 


-1.4 

-  2  I  max  .  f 

-  OTW  /  YACOB(I.K)  •  DELTAT(I.K) 

-  00(11,11. 1-1, K)  -  0T1  •  YAC0B(1-1  ,K) 

-  EE\ 1 1 ,11,1-1 ,K)  -  DTI  •  YAC08( 1+1 ,K) 

-  1.  +  2.  •  OTW  •  DELTAT(I.K) 

-  2  .  I MAX  -  1 


GG 

[1.I-1.K 

■ 

0Q1{ 

[I.K 

• 

OELTAT 

[I.K) 

GO 

2.I-1.K 

m 

002 

I.K, 

i  • 

OELTAT 

I.K 

GG 

3.1-1  ,K 

m 

003 

I  ,K 

I  • 

OELTAT 

I.K 

GG 

[4.I-1.K 

0041 

I  ,K 

• 

OELTAT 

I.K 

999  CONTINUE 
3  CONTINUE 


PERFORM  BLOCK  TRIO  I  AGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 

CALL  MATRXI(IMAX.KMAX) 

DO  Ml  K  -  2  .  KMAX  -  1 


DO  991 

-  2 

IM1 

001 

I.K, 

-  GG 

.  i-1  .k; 

DQ2 

I.K, 

-  00 1 

,2.1-1  .K 

DQ3 

I.K’ 

-  GG 

,3,1-1  ,K 

DQ4{ 

I.K 

-  GGI 

4.I-1.K, 

991  CONTINUE 

K-SNEEP  BEGINS  HERE. 


DO  13  I  -  2  .  IM1 

CALL  AMAT2(I  .KMAX-1  .ZETAX.2ETAZ.ZETAT) 

DO  IS  II  -1.4 

00  IS  12  -1.4 

DO  15  K  -  2  .  KMAX  -  1 

EE(I1 , 12, 1  ,K-1)  -  A(1 1 . I2.K+1 )*DTH  •  DELTAT(I.K) 
00(11 , 12. 1 ,K-1 )  -  — A( II . 12, K— 1 ) *OTH  •  DELTAT(I.K) 
IS  CONTINUE 

C 

C»»»  SECOND  ORDER  DAMPING  AOOEO  HERE. 

C 


C 

C 

C 


DO  16  II 
DO  16  K 
0T1 

00( 11,11,1, K— 1 
EEC  1 1 ,11,1 ,K-1 
16  M4(I1, 11. I.K-1 


00 

17  K 

-  2  , 

KMAX 

GG 

[1. I.K-1] 

-  001 

[I.K) 

GG 

2. I.K-1 

-  002 

I.K) 

GO 

3. I.K-1 

-  0031 

I.K 

GGI 

[4. I.K-1] 

-  0041 

[I.K) 

-1.4 

-  2  ,  KMAX  -  1 

-  OTW  /  YACOB(I.K)  •  OELTAT(I.K) 

-  00(11, II. I.K-1)  -  DTI  •  YACOS(I.K-I) 

-  EE(I1, II. I.K-1)  -  0T1  •  YACOB( I ,K+1 ) 

-  1.  ♦  2.  •  OTW  •  OELTAT(I.K) 

1 


17  CONTINUE 
13  CONTINUE 


PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 


CALL  MATRX2(IMAX,KMAX) 

00  18  I  -  2  ,  IMAX  -  1 
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DO  18  K 
001(1. K) 
002(1, X) 
003(1. K) 
004(1. X) 
18  CONTINUE 


■  2  ,  Kill 

*  wn.i.K-i 

■  2.I.K-1 
-  WI3.I.K-1, 

*  M(4.I.X-1) 


C«*«  UPDATE  FLOW  VAWA8t.ES  AT  INTERIOR  POINTS. 
967  CONTINUE 

■AM  .  8. 

RUtfJT  -  «. 

MU  .  «. 

-  6. 

A  -  2  ,  KM1 

I  •  2  .  1M1 

-  on 

-  02 

-  03  ( 

-  04  ( 


qo  ns 

00  19 
01 (I.X 
I  ,K 

03(1. « 


(I.K 

♦  001 

(I.K. 

•  YAC08 

[I.K) 

I.K, 

+  002 

I.X 

•  YACOB 

I.K) 

(I.K 

+  0031 

'I.X 

•  YACOB 

I.X) 

(I.K] 

♦  D04| 

I.X) 

•  YACOB I 

I.X) 

19  CONTINUE 

Cl  1  M  lQQT0^Iyg  '"HERE  iV“5J2“®  DENSITY  CHANGING  MOST  RAPIOLY 

IF  (RMAX. LT.A8S(001  ( i  , K)»YAC08(I  ,K)))  THEN 

IR  •  I 

KR  .  K 

END  IF 

"MAX  -  AMAX1  (RMAX,ABS(DQ1  (I  ,K)  .  YAG0B(I.K)1 

“  AMAXljR^X.ABS(DQ2(I.K)  .  YAC0B(i.xj 

I  »K)  • 


-  AMAX1 C RVMAX, ABS(DQ3( 


YACOB(I.K) 


-  AMAX1  (EMAX.A8S (004(1  ,K)  •  YAC06(I ,K) 


RUMAX 
RVMAX 
EMAX 

999  CONTINUE 
e  ""te 

CIIIIISELECT  INTERVAL  AT  WHICH  OUTPUT  OF  RESIDUALS  IS  DESIRED 

1EMAx!ir’kR/1#*1*'EQ  (,™"1))  mi€  "44X.8DMAX. RVMAX, 

RETURN  ' 

3##2r  IR^SX^KR^j  ,W,MAX,  • 1 1X<  ’OUMAX’ .  1 1X.  ’DVMAX’ .  1 1X.  ’DEMAX’ ,  18X. 

3961  F0RUAT(4(E14.8.2X) .213) 

ENO 

SUBROUTINE  MATRX1  ( IMAX.KMAX) 

iccW8ni?)00<4,4,161,41),l*<(4’4'iett4,)>EE(4’4’1#,’4,>’ 

C0MM6jVSCRAT/A(4.4.161).HH(4,4.181.41),C(4,5,181) 

2,ut.tJ::H!:tJ!l4,’U2U2L42’LMLWL44 


00  1  II  -  1.  4 

00  1  X  -  2  .  KMAX  -  1 


00(11. 1.K)  •  AI 
EE(I1.1.1.X 


(M.1.X) 

(11.1.1  .XI  __  . 

11,2,1 , XI  .  EEC  II ,2.1,K 
(II. 3.1.X)  -  EEC  11.3,1 . X 
(11.4.1 .X)  -  EE (II ,4, 1  ,K 
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1  CONTINUE 


00  1«M  I  -  2  .  I  MAX  -  2 

00  5  It-  1  .  4 

00  2  X  •  2  .  KMAX  -  t 

c(ii.i.k)  -  gc(ii.i.k)  -  oo 

1  -  DO 

2  -  00 

3 _  -  DO 

2  CONTINUE 
OO  5  12  -  1  ,  4 

00  5  K  -  2  .  KMAX  -  1 

^A(I1 , 12, K)-  144(11 , 12, 1  ,K)  - 

2  I 

3 

C( 1 1 , 12+1 ,K)»  EE(I1 , 12, 1  ,K) 

9  CONTINUE 

00  3  K  -  2  ,  KMAX  -  1 
L1 1  -  A(t,t.K) 


1 1.1. 1.  X)  *  iXt 

11.2.1 , X  i  00 

11.3.1 , X 

11.4.1. x 


■  00 
•  00 


1.1- 1.X 

2.1- 1.K 

3.1- 1.X 

4.1- 1.K 


00 

11.1.1.x] 

■  W 1.12. 1-1 .X) 

DO 

,11.2.1.x 

*  *i(2, 12, 1-1  ,x' 

00 

,11.3.1  ,K 

*  Hh(3.I2.I-1.X' 

001 

1 1 , 4 , 1 .  X] 

•  *(4. 12. 1-1. X) 

1 .  /  L22 


*  U12 

•  U13 

•  U14 


(4(2. 3.X)  - 

(4(2. 4. X)  -  _ 

L32  -  AI3.2.X)  -  '.31  •  U12 
L42  -  AI4.2.K J  -  L41  •  U12 
L33  -  A(3,3,K)  -  L31  •  U13  - 


!: 


L2I 

L2I 


L3I 
U34 
L43 
L44 
L4I 
C 
C 
C 


1.  /IM  -  ‘ 

(4(3, 4, K)  -  L31  •  U14  -  132  •  U24) 

4(4, 3, K)  -  L41  •  U13  -  142  •  U23 
*.4, K)  -  L41  •  U14  -  L42  •  U24  -  L43 


1.  /  L44 


•  L3I 

*  U34 


1 


(i.i. x)  •  c(i.i.x)  • 
'2.1.X>  -  (C(2.1.X)  - 
(3.1.X)  -  (C(3.1,X)  - 


C(4.1,K) 


LI  I 
L21 
L31 

,  ,  -  L32 

(C(4, 1 ,X)  -  L41 


1,1.2.x)  -  C(1 ,2,X)  . 

I  2.2.X)  -  (C(2,2.K)  - 
(3.2.X)  -  (C(3,2.K)  - 


C(4.2.X)  -  (C(4,2,K)  -  L41 
C(1 ,3,X)  -  C(1 ,3,X)  •  LI  I 


LI  I 
L21 
L31 
-  L32 


1 


•  Ci  I.I,*)}  •  L2I 

•  C  1.1,*) 

•  C  2,  ■  .X  j)  •  L3I 

■  «(»,»<)  -  L42  •  C(2.1.X) 

-  *  C(3, 1 ,X))  •  L4I 

•  C(1,2,xn  •  L2I 

•  Cl  1.2.x) 

•  C(2,2,x))  •  L3I 

•  C(1.2,X)  -  L42  *  C(2,2.K) 

-  L43  •  C(3,2,X))  .  L4I 


CI2.3.X)  -  C2.3.X)  -  L21  • 
C(3.3.X)  .  (C(3,3,K)  -  L31  * 
-  L32  * 


CiM.xS 


C(4,3.K)  •  (C(4,3,K)  -  L41  * 


(1.4.X)  -  C(1 ,4.X)  • 
(2.4.x)  •  (C(2.4.K)  - 
(3.4.X)  -  (C(3.4,K)  - 


LI  I 
L21 
L31 


-  -  „  L2I 

J.X 

1.3. X  )  •  L3I 

1.3. X  -  L42  •  C(2,3,K) 

-  LO  •  C(3,3,X))  «  L4I 


C(3 

[!:::!}’  •  u‘ 
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noo 


1 


C(4.4,K) 


Cl 

Ci 

c 


1 


C(4,3,K) 


(3.1.K) 

(2.1.K) 


C(M.K) 


(3.2.K) 

(2.2.K) 


C(1.2.K) 


(3.3.K) 

(2.3.K) 


C( 1 ,3,K) 


-  L32  •  C(2.4.K>)  •  L31 
(C(4.4.K)  -  141  •  C(1.4.K)  -  142  •  C(2.4,K) 

III  *  143  •  C(3.4.K))  •  L4I 

1.3. K1)  •  L2I 
t.S.K  i 

2.3.  K 

1.3. K 


:p.s.K)  -  cn.s.K)  • 

(2.3.Ki  -  fC(2,5,K)  - 
1(3, 5, K)  -  {cp.S.Kj  - 


L21  •  C 
131  •  C 
-  132  •  C 
(C(4.3.K)  -  L41  •  C 


C(3.1,K)  -  U34 
C(2, 1 ,K)  - 


:(3,2.K)  -  U34 
5(2. 2. K)  -  U24 


U24 

,  -  U23 

C(1.1.K)  -  U14 

-  U13 

S .  . .  . 

-  U23 
C(1 ,2,K)  -  U14 

-  U13 
0(3. 3, K)  -  U34 
C(2,3,K)  -  U24 

,  -  U23 

C(1,3,K)  -  U14 

-  U13 

-  U34 

-  U24 

-  U23 


)  «  L3! 

.  -  L42  •  C(2.5,K) 

.  ;  ^43  •  C(3,3,K))  .  L4I 


C(3,4,K) 

-  CI3.4.K) 

C(2,4,K) 

-  C(2,4,K; 

C(1.4.K) 

-  C(1.4.K) 

C(3. 5.K) 

-  CI3.3.K) 

C(2,S,K) 

| 

-  C(2.5.K) 

C(1.3.K) 

-  C(1.3.K) 

1 

3  CONTINUE 

00  6  II 
DO  9  K 
9  CG(II.I.K) 
DO  6  12 
DO  6  K 


-  U13 


•1.4 
•  2  .  KMAX  -  1 

-  C(II.I.K) 
■1.4 

-  2  .  KMAX  -  1 


-  U12  •  C(2, 1 ,K) 


-  U12  •  C(2,2,K) 


-  012  *  C(2,3,K) 


-  U12  •  C(2,4,K) 


-  012  •  C(2,5,K) 


iee«  continue 


backward  substitution 

00  7  I  .  IMAX  -  3. 

00  7  II  -1.4 

00  7  K  *  2  .  KMAX  -  1 
,00(11. 1.K)  -CC(H.I.K)  -  HH 


1.-1 


-  HHI 

[I1.1.I.KJ 

1  •  OGI 

1.I+1.K 

-  HH| 

I1.2.I.K 

1  •  0G| 

2.I+1.K 

-  HH| 

11. 3.1, K 

1  •  CC| 

3.I+1.K 

-  HH( 

'll. 4.1. K) 

1  •  GG( 

4.1+I.Kl 

7  CONTINUE 
RETURN 
ENO 

SUBROUTINE  MATRX2( IMAX, KMAX) 

1GC(40^,e1R41)D0^4,4,,61•4,),l*,(4,4•1#’’4,)’EE^4'4•,6,>4,). 
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oooo 


C0*«I/9CRAT/A<4,4,iei).HH<4.4.161.41).C(4.9,161) 

REAL  m 

REAL  L1 1 . L21 , LSI . L41 . L22 . LS2.L42 , LSS.L4S.L44 
2.L1I.L2I.LSI.L4I 

THIS  SUBROUTINE  PERFORMS  THE  BLOCK  TRIDI  AGONAL  MATRIX  I  VERSION  FOR 
AN  ENTIRE  ^CONSTANT  PLANE  DURING  THE  ZETA-  SWEEP 


DO 

1  11 

-  1. 

4 

DO 

1  I 

-  2 

IMAX  -  1 

AI 

-  1. 

/  I6fl(1,1, 

OG 

;n .  i.D 

-  GG 

11. 1.1)  • 

HH 

ii. i.i.i] 

-  EE 

n.  i.i.i; 

HH 

11.2.1.1 

-  EE 

,11.2.1.1 

HH 

11.3.1,1] 

-  EE 

11 .3.1.1 

HH! 

[11.4,1.1) 

-  EE 

11.4,1.1 

1  CONTINUE 


AI 

•  AI 

•  AI 

•  AI 

•  AI 


DO  IBM  K  •  2  .  KMAX  -  2 

DO  S  II-  1  .  4 

DO  2  I  -  2  .  IMAX  -  1 

C(II.I.I)  -  CC(II.I.K)  -  DO 

1  -  DO 

2  -  DO 

3  -  DO 

2  CONTINUE 

DO  S  12  -  1  .  4 

DO  S  I  -  2  .  IMAX  -  1 

A(I1,I2.I)-  161(11. 12. 1. K)  -  00 
1  -  00 

2  -DO 

C( I 1,1241. I)-  EE(I1 , 12, I ,K) 

9  CONTINUE 


I1.1.I.KI  ■  MO. I. K-1 

I  4  A  •  U  u  J  A  I  U  J 


11 .2.1 .  K 

11.3.1 . K 

11.4. 1. K 


*  GO 

*  00 
*  00 


2.1. K-1 

3.1. K-1 

4.1. K-1 


11. 1. 1. K)  ■  mil, 12,1. K-1 

11 .2.1. Ki  *  mt2.I2.I.K-1 

11.3.1. KS  i  mfj.  12, 1 ,K-1 
1 1 ,4, 1  ,K)  *  mi|4.I2.I.K-1 


00  3  I  -  2 
L1 1  -  A(l.l.l) 
LI  I  -  1.  /  L1 1 


IMAX  -  1 


LI  I 
LI  I 
LI  I 


-  L21  •  U12 


L2I  -  1 .  /  L22 


U23  -  (A(2.3.I)  -  L21  •  U13) 
U24  -  (A(2,4, I)  -  L21  •  U14) 
L32  -  A(3,2, I)  -  L31  •  U12 
L42  •  A(4.2.I)  -  L41  •  U12 


L2I 

L2I 


L33  -  A(3.3.I)  -  L31  •  U13  -  L32  •  U23 
L3I  -  1 .  /  L33 

U34  -  (A(3,4, 1 )  -  LSI  •  U14  -  L32  •  U24)  •  LSI 

L43  -  A(4,3, l)  -  L41  •  U13  -  L42  •  U23 

L44  -  A(4,4,l)  -  L41  •  U14  -  L42  •  U24  -  L43  •  U34 


/  L44 

LI  I 

L21  •  C 
L31  •  C 

1  -  L32  •  C 

C(4, 1,1)  -  (C(4, 1,1)  -  L41  .  C 


L4I  -  1. 
C 
C 
C 


1.1.1)  -  C(1.1,I)  .  I 

2.1.1)  -  (C(2,1.I)  - 

3.1.1)  -  (C(3, 1,1)  - 


1.1.1) )  •  L2I 

1.1.1) 

2.1.1) )  •  LSI 

1.1.1)  -  L42  •  C(2,1 , I) 
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-  143  *  C(3, 1.1))  *  L4I 


c  1,2.1 

C  2,2,1 
0(3.2, I 

- 

- 

- 

0(1, 2.1)  »  LI I 
(C(2 , 2,1)  -  L21 
(0(3,2, I)  ~  L31 

1 

-  L32 

0(4.2. I)  - 

(C(4.2. I)  -  L41 

C(  1 ,3. 1 

- 

C(1 ,3, 1)  *  HI 

0(2.3. I 
0(3,3. I 

- 

- 

(C(2 , 3,1)  ~  L21 
(0(3 , 3,1)  ~  L31 

-  L32 

C(4,3, I )  - 

(C(4,3, I)  -  L41 

0(1.4.11 
0(2.4, I 
0(3.4, I) 

- 

C(1 ,4, 1)  •  LI  I 
(C(2, 4,1)  -  L21 
(0(3, 4,1)  -  L31 

-  L32 

0(4,4. I)  - 

(0(4.4, I)  -  L41 

0(1, 5, I) 
0(2.5. I 

0(3.5. I) 

* 

• 

0(1. 5.1)  *  LI  I 
(0(2.5, I)  ~  L21 
(0(3,5, I)  -  L31 

-  L32 

0(4,5. I)  - 

(0(4.5. I)  ~  L41 

0(3, 1.1) 

m 

0(3, 1.1)  -  U34  * 

0(2, l.l) 

■ 

0(2.1  ,l)  -  U24  * 

-  U23  * 

0(1. 1.1) 

■ 

C( 1 . 1  . 1)  -  U14  • 

-  U13  * 

0(3.2. I) 

■ 

C(3,2, I)  -  U34  • 

0(2,2. 1) 

* 

0(2,2. l)  -  U24  « 

C(  1,2.1) 

1 

C(3,3, I ) 
C(2.3. I) 

1 

C( 1 ,3, 1 ) 


C(3,4, I ) 
C(2,4. I) 

1 

C(  1 1 4 , I ) 

1 

0(3.5. 1) 
0(2.5. i) 


[1.2. 1) )  •  L2I 
1.2.1) 

[2.2.1) )  •  L3I 

[1.2. 1)  -  L42 

-  L43  •  C(3, 

[1.3.1) )  *  L2I 

1.3.1) 

2.3. 1) ) 

[1.3. 1)  - 

-  143  * 


L3I 

L42 

0(3, 


5(1.4,!))  •  L2I 
:( 1 .4.  ] 


•  C 

»  C( 1 ,4. I ) 

•  C  ( 2 , 4 , 1 )  )  •  L3I 

•  C(1 ,4, I )  -  L42 

-  L43  •  C(3, ■ 


C(3,3, I )  - 
C(2 . 3 . I )  - 


-  U23 
C(1 .2. I)  -  U14 

U13 
U34 
U24 
U23 

C(1 .3. I)  -  U14 

-  U13 
C(3.4, I)  -  U34 
C(2.4. I)  -  U24 

-  U23 
C( 1 .4.1)  -  U14 

-  U13 
0(3,5, I)  -  U34 
C(2 ,5,1)  -  U24 

-  U23 

0(1. 5,1)  -  C(1 ,5, I)  -  U14 

U13 


3  CONTINUE 


C(1 ,5, I))  •  L2I 
C  1.5.1) 

C(2,5, I))  .  L3I 
C(1 ,5. I)  -  L42  ' 
-  L43  •  C ( 3 . ! 
4.1.J) 


-  U12  • 


3.2.1) 

4.2.1) 

3.2.1)  -  U12  * 

4.3.1) 

4.3.1) 

3.3.1) 

4.3.1) 

3.3. 1)  -  U12  • 

4.4.1) 

4.4.1 
3,4,0 
4.4,0 

3,4.0 

4.5.1) 

4,5,0 

3.5.1 

4.5.1 

3.5. 1 


-  U12  • 


-  U12  • 


DO  6  II 
DO  9  I 

9  GG(II.I.K) 
DO  6  12 
00  0  I 
HH(I1 ,12.1 
6  CONTINUE 
1000  CONTINUE 


-1.4 

-  2  ,  I  MAX  -  t 

-  C(ll.l.l) 

-1.4 

-  2  ,  IMAX  -  1 
K)  -  C(I1. 12+1,1) 


•  C(2,2, I) 

2,0)  •  t-4I 

«  C(2.3.I) 
3.1))  •  L4I 

•  C(2,4, 1) 
4.0)  •  L4I 

•  C(2 . 5, I ) 
5.0)  •  L4I 

C(2. 1,1) 

0(2,2, I) 

C(2,3, I) 

C(2,4, I) 

C(2,5, I) 
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BACKWARD  SUBSTITUTION 


DO  7  K  -  KMAX  -  3,  1  .  - 

1 

DO  7  11  -1.4 

DO  7  I  -  2  .  IMAX  -  1 

GG(ll.I.K)  -  GG(II.I.K)  -  HH 

1 1 1 , 1 , 1 ,K , 

*  GG 

(1.I.K+1) 

1  -  HH 

I1.2.I.K 

.  GG 

2.I.K+1) 

2  -  HH 

I 1 ,3, I ,K] 

*  GG 

3.I.K+1) 

3  -  HH| 

11.4,1 ,K] 

*  GG( 

4.I.K+1) 

7  CONTINUE 


RETURN 

ENO 

SUBROUTINE  METRIC 
COMkION/F  I X/OMEGA 

COMMON/DGRID/DT, IMAX.KMAX,  ITEL.ITEU 
C0fcM0N/GRID1/X(161 , 41 ) ,2(161 ,41) 

COMkON/GR I O/YACOB (181 ,41) 

COMAON/MTR I X/X i X ( 161 ,41),XIZ(161,41),  ZETAX(  161 .41 ) . ZETAZ( 161 ,41 )  , 
1XIT(161 ,41), ZETAT (161,41) 

••  SUBROUTINE  METRIC  COMPUTES  THE  METRICS  IN  BOTH  DIRECTIONS  AND 
THE  UNSTEADY  COEFFICIENTS  ETAT,  ETC. 


DO  1000  K  -  1 
DO  1000  I  -  1 
XTAU  -  OMEGA  • 

YTAU  -  OMEGA  • 

C**»  PRESENT  SET  UP  IS  FOR  FI.OW  PAST  AN  AIRFOIL. 

C 

C! ! I  I (CENTRAL  DIFFERENCES  AT  INTERIOR  POINTS,  TWO-POINT  ONE-SIDED 


,  KMAX 
.  I  MAX 
Z(I,K) 

(-X(I.K)) 

IS  FOR  FLOW  PAST  AN 


C! 


—  — -  -  •  •  — •  '  -  — >  ■  »  *  ■  ■  ■  W  I  ■  I  w  I  IT  r  VMk  —  1  VLU 

INDIFFERENCES  DOWNSTREAM,  THREE-POINT  AT  OTHER  OUTER  BOUNDARIES 
IF(I.EQ.I.OR.I.EQ.IMAX)  GO  TO  10 
XXI  •  .5 
ZXI  -  .5 
GO  TO  15 

10  IF(I.EQ.IMAX)  GO  TO  16 

XXI  -1.0.  (X(2 . K )  -  X(I.K)) 

:(2.K)  -  Z(1. 


l  .  un  .  1  .  CU  .  1 MAA  J  W  IU  1 

*  (X( 1+1 . K ) — X (1—1 ,K)) 

.  (Z(I+1,K)-Z(I-1,K)) 


(Z( 


K)) 


17 


ZXI  -  1.0 
GO  TO  15 

16  XXI  -  1.0  .  (X( IMAX , K)  -  X( I MAX 
ZXI  -  1.0 
15  CONTINUE 
IF(K.EQ.1 
XZET  -  .5 
ZZET  -  .5 
GO  TO  20 

IF(K.EQ.KMAX)  GO  TO  18 
XZET  -  2.  *  X' 


(Z(IMAX.K)  -  Z( IMAX 

. OR.K.EQ.KMAX)  GO  TO  17 
»(X(I,K+1)-X(I,K-1)) 
»(Z(I.K+1)-Z(I,K-1)) 


K:!}{ 


ZZET  -  2. 

GO  TO  20 

18  XZET  -1.5.X 
ZZET  -  1 .5  »  Z 
20  CONTINUE 
Cl  I  I  I ICOMPUTE  JACOBIAN 

YACOBI  -  XXI  .  ZZET  -  XZET  . 
YACOB(I.K)  -  1.  /  YACOBI 
XIXfl.K)  -  ZZET  •  YACOB(I.K) 
XIZ(I.K)  -  -XZET  •  YACOB(I.K) 
XTAU  -  OMEGA  .  Z(I ,K) 


(1.2) -15  .  X ( I , 1 )  -  .5  .  X( I .3) 

(1.2)  -  1.5  .  Z(I,1)  -  .5  .  Z(I ,3) 

( I , KMAXI-2 . •  Xf  I  .KMAX-1 )+.5«X( I .KMAX-2) 
(I ,KMAX)-2.  *  Z( I .KMAX-1 )+. 5.Z( I . KMAX-2) 


ZXI 
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*  -  '  nvvu  \  i  ,  r\  f 

ZETAX(I.K)  .  XTAU  -  ZETAZ(I.K)  .  YTAU 


YTAU  -  -  OMEGA  .  X(l,K) 

"  "  XIX(l.K)  .  XTAU  -  XIZ(I.K)  .  YTAU 
ZETAX(l.K)  -  -ZXI  •  YACOe(I.K) 

ZETAZl I ,K)  -  XXI  .  YAC06( I  ,K) 

ZETAT(I , K )  -  -  -  -  ■ 

1000  CONTINUE 
RETURN 
END 

SUBROUTINE  DISSIP 

COfAION/FLOW/QI  {161 ,41), 02 (161 , 41 ) ,QJ( 181 ,41)  04(161  411 
COMMON/PERTR/DQ1 (161,41)  .002(161 .41 ), 003(1 61  41 )  004(161  41) 
CO^WON/DGRID/OT.IMAX.KMAX.ITEL.ITEU  ^  U  } 

CO0ON/GR I D/YACOB  (161,41) 

COfcWON/D AMP/WW , WW I 

DIMENSION  P(161),EPS(161),DIS1(161 . 4) ,DIS2( 1 61 .4) 

SfiSrlSKS  AD0S  ™E  F0URTH  0R0ER  DISSIPATION  TERMS  TO  THE 

IM1  -  IMAX  -  1 
KM1  -  KMAX  -  1 
I M2  -  IMAX  -  2 
KM2  -  KMAX  -  2 

DO  10  K-2  ,  KM1 

COMPUTE  SWTICHING  FUNCTION  BASED  ON  SECOND  DERIVATIVE  OF  PRESSURE 
w  1  I  ■  1  ,  I MAX 

1  00*2  I-?  *  jJJJI-K)-(Q2(,.K)‘*24Q3(I,K)..2)/(2..Q1(I.K))) 

IP2  -  I  +  2 

IF(I .EQ. IM1 )  IP2  -  IMAX 
IM2  -  I  -  2 


I F ( I . EQ . 2 )  I M2  -  1 
IP1  -  I  +  1 
IF(I.EQ. IMAX)  IPt  . 
IM  -  I  -  1 
IF(I.EO.I)  IM  -  1 
IF(I.EQ.I)  I M2  -  1 
IF(I.EQ.IMAX)  IP2  - 


IMAX 


IMAX 


x:’?"  !•«>«*&<  i;?!#r<p(  ,p'  >«•  -p<  ■  »«■(  ■“» 

DISin.i)  -  (Q1  (IP1  ,K)-Q1  ( I  ,K))»VOL 

°  !  ){-2j  “  (Q2( IP1 ,K)-Q2( I ,K))«VOL 

DISUI.J)  -  (Q3( IP1  , K)-Q3(  I  ,K))»VOL 

HP1  -  04  ( IP1 ,K)+P(IP1) 

HP  -  Q4( I ,K)+P( I ) 

HM1  -  Q4IIM.K)  +  P(IM) 

HP2  -  Q4( IP2.K)  +  P( IP2) 

HPM  -  04(IM,K)+P(IM) 

0IS1 ( I , 4)  -  (HPI-HP).VOL 

lih'i  r  te3i;:js}K:!®}!:*{fSjiS:SJsa 

Sia  !:J  - 

CONTINUE 

DO  15  I  -  1  .  IM1 

02P  -  AMAX 1 ( EPS( I ) . EPS( 1+1 ) ) 

C22  -  60.  .  02P 

C11  -  AMAX1 (0. 0, ( 1  -C22)) 
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C22  -  C22  •  WWAACOB(  I  ,K)  «  OT 
C11  -  C11  •  WWAAC06(I  ,K)  •  DT 

CIIIHSWITCH  ON  SECOND-ORDER  AND  SWITCH  OFF  FOURT:4-ORDER  DISSIPATION 
CHIMIN  VICINITY  OF  SHOCKS 

DISI(I.I)  -  C11  «  DIS2(  1,1)  +  C22  •  DISI(I.I) 

DIS1  (1 ,2)  -  C11  *  DIS2( I , 2 )  +  C22  *  01SUI.2) 

OIS1 (1.3)  -  C11  •  DIS2( I ,3)  +  C22  •  DIS1CI.3) 

DIS1 (1,4)  -  C11  •  DIS2( I ,4)  +  C22  •  0ISl(l,4) 

15  CONTINUE 

DO  16  I  -  2  .  IM1 


D01( 

I.K 

-  DIS1 

1,1) 

-  DIS1 

1-1.1) 

DQ2 

I.K 

-  0IS1 

1,2' 

-  DIS1 

1-1,2 

003 

I.K 

-  DIS1 

1,3 

-  DIS1 

1-1.3) 

DQ4{ 

I.K 

-  DIS1 

1.4) 

-  D I S 1 

1-1,4) 

16  CONTINUE 
10  CONTINUE 


K  DIRECTION 

!  I  I  I ! FOURTH-ORDER  DISSIPATION  ONLY 
DO  30  I  -  2  .  IM1 
WT-  0.5  •  DT  •  WW  /  YACOB( I ,2) 

W3  -  0.5  •  DT  •  WW  /  YACOB(I.KMI) 

DQ1 (1,2)  -WT.  (Q1(I,1)  -  2.  •  01(1,2)  +  Q1(I ,3)) 
1+001(1,2) 

002(1,2)  -WT.  (02(1,1)  -  2.  »  02(1,2)  +  02(1.3)) 
1+002(1,2) 

DQ3(I ,2)  -WT.  (03(1,1)  -  2.  .  03(1,2)  +  03(1,3)) 
1+003(1,2) 

DQ4(I ,2)  -WT.  (04(1.1)  -  2.  .  04(1,2)  +  04(1.3)) 
1+004(1.2) 

WT-  W3 

DQI(I.KMI)  -WT.  (Q1 ( I ,KM2)  -  2.  .  01(1, KM1)  +  01 ( I ,KMAX) ) 
1+001(1, KM1) 

002(1, KM1)  -WT.  (Q2(I,KM2)  -  2.  *  02(1, KM1)  +  02(1, KMAX)) 
1+DQ2( I ,KM1 ) 

003(1, KM1)  -WT.  (03(1, KM2)  -  2.  *  Q3(I.KM1)  +  Q3(I,KMAX)) 
1+0Q3(l,KM1) 

DQ4( I ,KM1 )  -WT.  (04(1, KM2)  -  2.  •  Q4(I,KM1)  +  04(1, KMAX)) 
1+DQ4(I.KM1) 

DO  35  K  -  3  .  KM2 

WT-  -  DT  .  WW  /  YACOB(I.K) 

DQI(I.K)  -WT.  (Q1 ( I ,K+2)  -  4.  .  01(1, K+1)  +  6.  *  Q1 ( 

II. K)  -  4.  .  Q1 ( I ,K-1 )  +  Q1(I.K-2))+DQ1(I.K) 

DQ2 ( I , K )  — WT.(Q2( I ,K+2)  -  4.  .  02(1, K+1)  +  8.  *  Q2( 

1I.K)  -  4.  .  Q2( I . K— 1 )  +  02(1  ,K-2) )+DQ2( I ,K) 

DQ3(  I  ,K)  -WT.(Q3(I  ,K+2)  -  4.  .  Q3(I.K+1)  +  «.  *  03( 

II. K)  -  4.  .  Q3(I.K-1)  +  03  ( I ,  K— 2 ) )+0Q3( I ,K) 

0Q4(I,K)  -WT.(Q4(I,K+2)  -  4.  •  04(1, K+1)  +  6.  *  Q4( 

11. K)  -  4.  .  04(1, K-1)  +  Q4(  I ,K-2) )+0Q4( I ,K) 

35  CONTINUE 
30  CONTINUE 
C 

RETURN 

END 

SUBROUTINE  WALL8C 
C0M40N/SURF/PSUR( 161) 

C0U»N/GRID1/X(161  .41 )  .2(161 .41) 

C0M40N/PAR/GAMMA ,  REYREF  .ALFA ,  ALFA1  .REDFRE,  AMINF ,  ALFAI 
COMMON/DGRID/DT , IMAX , KMAX , I  TEL, ITEU 
C0M40N/GR I D A ACOB ( 1 6 1 . 4 1 ) 

C0M40N/DAMP/WW ,  WW  I 


98 


Cl  1 1 


.St¥oS<YI?)*iE^(?4;*I?ix,z<’**-4’)zeTAX<’*,*,)-2ETAZ(,e,-4')- 

41  }  02O61  .♦D.OJde1 .4’ ) .04(161  .-»1 ) 

u I MENS I ON  C 1 ( 4 ) 

DIMENSION  A(2.2).RHS(2) 

I  I  APPLY  BOUNDARY  CONDITIONS  ON  THE  CUT  AND  THE  AIRFOIL  SURFACE 
DO  9  I-ITEU, IMAX 


11 

■ 

IMAX  4 

1  - 

Q1 

(I. 

1) 

■  . 

5  • 

02 

I. 

1) 

-.5 

•  (02 

03 

I. 

1) 

m 

5  • 

04 

I. 

1) 

-  .5  . 

01( 

*11 

.1 

1-01 

II. 1 

02 

11 

.1 

1-02 

*1.1 

03 

11 

.1 

)-Q3 

1.1 

04 

11 

.1 

■04  ( 

1.1 

(01(1.2)401(11,2)) 
(1,2)402(11,2)) 
(03(1,2)  4  03(11,2] 
(04(1,2)404(11.2)) 


ITEU 


UC0N3  -  (02(l ,K)«C1 (2)403( I ,K)»C1 (3) ) 
1/Q1 ( I ,K) 

K  -  2 

C1(1)  -  XIT(I.K) 

Cl (2)  -  XIX(I.K) 

Cl (3)  -  XIZ(I ,K 

1/QUI  -K)*C1 (2)403(1  ,K)»C1(3)) 

RHS(1)  -  2.  •  UC0N2  -  UC0N3  -  XIT( 1.1) 
FOR  VISCOUS  FLOWS  SET  UCON  TO  ZERO  ALSO 
IF(REYREF.GT.0. )  RHS( 1 )  -  -  XIT( I , 1 ) 

*'*  *'  XIXjI.1)  ' 

-  ZETAX(I.I) 

-  ZETAZ(I.I) 

.  -  -  -  ZETAT (1,1) 

TEW3'  -  A ( 1  , 1 ) 

IEW2  -  A(  1 , 2) 

T&FJ  *  A(2, 1 ) 

TQ*i  m  a(2.2) 

DEN  *  1 


*(*.*)  * 

-  XIZ(l 


, ,  /(TEMPI  •  TEMP4  -  TEMP 2  .  TEMP3) 

*l.(  -  A(2, 2)  «  DEN 
*0  1)  -  -  TEMP2  •  DEN 
AU.t I  *  -  TEMP 3  .  DEN 
A(a,2l  m  TEMPI  •  den 
01(1.1)  -  2.  «  01(1,2)  -  01(1,3) 

0.1  ■  01(1.1  •(A(1,1).RHS(1)4A(1,2).RHS(2)) 

1  CONtInue"  01  *(M2 , 1 ) »RHS(  1 )4A(2 , 2) »RHS(2) ) 

DO  10  I-ITEL  , ITEU 
U2-Q2( 1 1 2)/Q1 (1,2) 

W2-03(I,2)/01(I,2) 

U>Q2(1L3)/O1(IQ3)I'2)"0  5*Q1  ( 1 .2)»(U2.U2+W2.W2)) 
W3-Q3(l!3)/01(li3) 

W04*  1  •  3)-0 . 5*01  ( 1 . 3)  »(U3»U34W3»W3)  ) 
P1**(4. «P2-P3)/3. 

PSUR(  I  )-(GA*41A.P1-1  .  )/(  .  7.AMINF..2) 
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U1-02fI.il/Q1 (1.1) 

W1-Q3(I.1)/Q1(I.1) 

16  04(1,1  )-P1 /(GAMMA— 1 .  )+0 . 5*Q1  (l,1)»(U1  *U1+#1  »*»1 ) 

RETURN 

ENO 

SUBROUTINE  STRESS ( I TN.OALFA) 

COJ0CN/FLO0/Q1  (161, 41),  Q2(161. 41)  ,03(161 ,41 )  ,Q4(  161  .41) 
COWON/DGRID/OT,  IMAX.KMAX,  ITEL.ITEU  ' 

COMMON/GR I D 1 /X ( 161,41), Z(161,41) 

’ REYREF  • A<-FA .  ALFA1 ,  REDFRE .  AMINF .  ALFA  I 

DIMENSION  AA( 1 61 .41) 

1 ,RH1 (161 ) ,RH2( 161) ,RH3( 161 ) ,RH4( 161 ) 

COMMON/LOGIC/RSTRT , PITCH , PLUNGE 
LOGICAL  RSTRT, PITCH, PLUNGE 
U(I.J)  -  02(1, J)  /  01(1, J) 

V(I.J)  -  03(1. J)  /  Q1(|, j) 

SUBROUTINE  ADDS  VISCOUS  TERMS  TO  THE  RIGHT  HAND  SIDE 
GOGM  -  GAMA  /  (GAMMA  -  1  .  ) 

IF(ITN.GT.10.OR. (RSTRT) )  CALL  EDDY(DALFA) 

COMPUTE  U  AND  V 
KMAXM1  -  KMAX  -  1 
IMAXM1  -  IMAX  -  1 
PR  -  1  . 

DO  10  K  -  1  .  KMAX 
DO  10  I  -  1  .  IMAX 

„  ( 0,<’’K)  - • 5  •  w(i.k)**»»(..k)..2) 

COMPUTE  TXX.TXY  AND  VISCOUS  DISSIPATION  AT  I  -  1  /  2 


DO  30 
DO  20 
UXI  - 
VXI  - 
AX  I  - 
UZET- 
V2ET- 
AZET- 
XXI  - 
ZXI  - 
XZET- 
ZZET- 
YAC  • 
YAC  - 
XIX  - 
ZETAX- 
XIZ  - 
ZETAZ- 
CNM  - 
UX  - 


K  -  2  ,  KMAXM1 
I  -  2  .  IMAX 
U(I.K)  -  U(I-I.K) 

V(I.K)  -  V(I-1 ,K) 

AA(I,K)  -  AA( 1-1 ,K) 

'  I ,K+1  1  ,K~1  >+U( 1-1  ■ K+1 I"1  .K-1  )  ) 

•  25. ( V(  , K+1 )-V( I . K-1 )+V(  1-1 , K+1 )-v( 1-1 , K-1 ) ) 

x(  I  *  1  ,K_1  )+AA(  1-1  •K+1  )-AA( 1-1  .*<-1 )  ) 

zu!k)  -  z(i-i ;K) 

«  *  ,K+1)-X(I  .K-ll+Xfl-l , K+1 )— X(  1—1 . K— 1 ) ) 

1 .  /YAC 
ZZET  .  YAC 
•  -  ZXI  •  YAC 
-XZET  •  YAC 
■  XXI  *  YAC 

.5  *  (CMU(I.K)  +  CMU(I-I.K)) 

UXI  .  XIX  +  UZET  .  ZETAX 

VXI  *  XIX  +  VZET  •  ZETAX 

AX  I  .  XIX  +  AZET  .  ZETAX 

UXI  »  XIZ  +  UZET  •  ZETAZ 

VXI  •  XIZ  +  VZET  •  ZETAZ 

AX I  »  XIZ  +  AZET  •  ZETAZ 

-(-4.  .  UX  +  2.  *  VZ)  *  CNM  /  3. 

CNM  .  (UZ  +  VX) 

-CNM  /  3.  .  (-4.  .  VZ  +  2.  »  UX) 
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R4 


S4 


-  ( (U( I .K)+U( 1-1 , K) )*TXX+(V( I ,K)+V( 1-1 ,K))»TXY)»0. 5 
+  CNM  /  PR/(GAM4A  -  1 . )  .  AX 

-  ((U(I,K)+U(I-1.K)).TXY+(V(I,K)+V(I-1.K))»TYY)*0.5 
+  CNM  /  PR  /  (GAK44A  -  1 . )  »  AZ 

DEBUG 

TURN  OFF  ENRGY  DISSIPATION  AND  DIFFUSION 
R4  -  e. 


ii :  fS!!S 


RH2 
RH3 
20  RH4 

DO  30  I  - 
DOT  ( I 
002  (  I 
DQ3 
004 

30  CONTINUE 

IN  THE  2  DIRECTION 


•  TXX  +  XIZ 

•  TXY  +  XIZ 

•  R4  +  XIZ  • 
IMAXM1 


)  /  YAC 
)  /  YAC 


•  TXY 

•  TYY 
S4)  /  YAC 


-  DQIfl.K)  +  RH1 1 

(I+’i 

\  -  RHl 

(i) 

-  D02( I ,K)  +  RH2I 

1+1 

\  -  RH2I 

] 

-  DQ3 ( I . K )  +  RH3I 

:i+ii 

1  -  RH3 

0 

-  DQ4( I ,K)  +  RH4( 1+1 j 

I  -  JW4| 

[[) 

I  - 


.25 

.25 

.25 

.25 


) 


DO  70 
DO  60  K  - 
UXI  -  .25 
VXI  - 
AX  I  - 
XXI  - 
ZXI  - 
UZET 
VZET 

AZET  -  AA(I.K)  -  AA(I.K-I) 
XZET  - 
ZZET 
YAC 
YAC  - 
XIX  - 
ZETAX 
XIZ  - 
ZETAZ 
CNM  - 
UX  - 


IMAXM1 

KMAX 

U ( I  + 1 . K ) — U (I  —  1 ,  K  )+U(  I  + 1 , K — 1 ) — U( I — 1 . K — 1 ) ) 

y(l+1  .K)-V( 1—1 ,K)+V( 1+1 , K-1 )-V( 1-1 , K- 1 ) ) 


AA(I  +  1  ,  K)-AA( 1-1 , K)+AA( 1  +  1 , K-1 )-AA( 1-1 .K-1)) 

X(I+1.K)-X(I-1.r .  . 

Z(I+1 ,K)-Z( 1-1 ,) 

:-D 


xn+i,K)-xn-i.K)+xn+i.K-n-x(i-i,K-i)) 

.  K)+Z(I+1.K-1)-Z(I-1.K-1)) 


-  U(I.K)  -  UfI.K-1 

-  V(I.K)  -  V( I , K-1 

-  AA( I ,K)  -  AA( I ,K 

-  X( I ,K)  -  X(I.K-I) 

-  Z(I.K)  -  Z(I.K-I) 

■  XXI  •  ZZET  -  ZXI  .  XZET 
|  1.  /  YAC 

ZZET  .  YAC 
'  -  ZXI  •  YAC 
-XZET  .  YAC 
XXI  .  YAC 
.5  .  (CMU(I.K)  + 


VX 

AX 

UZ 

vz 

AZ 

TXX 


-  UXI 

-  VXI 

-  AX  I 

-  UXI 

-  VXI 

-  AX  I 
(-4. 


TXY  -  CNM  • 


XIX  +  UZET 
XIX  +  VZET 
XIX  +  AZET 
XIZ  +  UZET 
XIZ  +  VZET 
XIZ  +  AZET 
•  UX  +  2.  • 
(UZ  +  VX) 


CMU(1 .K-1)) 

•  ZETAX 

•  ZETAX 

•  ZETAX 

•  ZETAZ 

•  ZETAZ 

•  ZETAZ 

VZ)  •  CNM  /  3. 


TYY  -  -CNM  /  3.  .  (-4.  .  VZ  +  2.  .  UX) 


R4 


S4 


R4  -  0 
S4  -  0 


((U(IK)+U(I.K-1)).TXX+(V(I,K)+V(I,K-1)).TXY).0.5 
+  CNM  /  PR/fGAAWA  -  1 . )  .  AX 
-  ((U(I,K)+U(I,K-1) ) »TXY+(  V(  I  ,K)+V( I ,K-1 ) )*TYY) *0 . 5 
+  CNM  /  PR  /  (GAMMA  -  1 . )  •  AZ  ’ 


RH1 


RH2(K) 
RH3 


K)  -  0. 

(ZETAX 
(ZETAX 
(ZETAX 


•  TXY)  /  YAC 

•  TYY)  /  YAC 


TXX  +  ZETAZ 
,  .  TXY  +  ZETAZ 

60  RH4(K)  -  (ZETAX  «  R4  +  ZETAZ  *  S4)  /  YAC 
DO  70  K  -  2  .  KMAXM1 

DQI(I.K)  -  DQI(I.K)  +  RH1 (K+1 )  -  RH1 (K) 
002 ( I , X )  -  0Q2 ( I , K )  +  RH2(K+1 )  -  RH2(K) 
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003 { 

I.K) 

-  003  ( 

I.K) 

+  RH3 

K+1 ) 

-  RH3(K) 

004  ( 

I.K) 

-  DQ4( 

I.K) 

+  RH4 

K+1 ) 

-  RH4(K) 

70  CONTINUE 
C 

RETURN 
END 

SUBROUT  I NE  LOAD ( CPS, CL. CD, CM. AL FAS) 

C0)nM0N/GRID1/X(  181. 4-1 )  .  Y(  1 61  .41) 

C0M4ON/DGRID/DT,  IMAX.KMAX,  ITEL.  ITEU 
DIMENSION  CPS( 101 ) 

THIS  SUBROUTINE  COMPUTES  THE  INVISCID  CONTRIBUTIONS 
TO  LOADS  ON  THE  AIRFOIL  SURFACE 

IMAXM1  -  IMAX  -  1 
CL  -  0. 

CO  -  0. 

CM  ■  0. 

DO  400  I  -  ITEL  .  ITEU  -  1 

xl  -  .5  •  (xn.n+xn+i.i)) 

YL  -  .5  •  (Y(1 , 1 )+Y( 1  +  1 ,  1 ) ) 

OX  -  X( 1+1,1)  -  X(I,1) 

DY  -  Y( 1+1,1)  -  Y(I.I) 

CPA  -  CPS(I+1)  •  .5  +  CPS(I)  *  .5 
DCL  -  CPA  «  (-OX) 

DCD  -  CPA  •  DY 
CL  -  CL  +  DCL 
CD  -  CD  +  DCD 

400  CM  -  CM  +  DCD  «  YL  -  DCL  •  XL 

DCL  -  CL  •  COS(ALFAS)  -  CD  •  SIN(ALFAS) 

CD  -  CL  «  SIN(ALFAS)  +  CO  •  COS(ALFAS) 

CL  -  DCL 
RETURN 
ENO 

SUBROUTINE  WRAP(  1 1 ,  J J  .  XSING,  YSING, XP,  YP.S0. AO.  Y0) 

DIMENSION  S0(161.4),Y0(41,4) , A0(  1 61 ,4) ,XP( 1 ) , YP( 1 ) 

THIS  SUBROUTINE  UNWRAPS  THE  AIRFOIL  AND  STORES  THE  UNWRAPPED 
SURFACE  IN  ARRAYS  A0  AND  SO.  IT  ALSO  DETERMINES  THE  STRETCHING 
IN  THE  ETA  DIRECTION. 

IMID  -  (II  +  1)  /  2 
DY  -  .8  /  (JJ  -  2) 

DO  1  J  -  2  ,  JJ 
Y  -  FLOAT (J-2)  •  DY 
1  YOfJ.1)  -  1 .25  •  Y  /  (1  .  -  Y.Y) 


Y0(  1  . 

D 

-  -  Y0(3, 1 

PI  -  4 

• 

ATAN  (  1  . ) 

ANGL  - 

PI 

+  PI 

U  -  XP| 

-  XSING 

V  -  YP| 

!d 

-  YSING 

U  -  1  . 

V  -  0. 

IIM1  - 

ii 

-  1 

DO  2  I 

m  1 

1  .  II 

XII  -  XP(I)  -  XSING 
Y11  -  YP( I )  -  YSING 

ANGL  -  ANGL  +  ATAN2{ (U»Y1 1-V«X1 1 ) ,  (U»X1 1+V»Y1 1 )) 
R  -  SQRT(X1 1 »»2  +  Y1 1 »*2) 
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U  -  XII 

V  -Y11 

D  •  CQOT /p\ 

A0(I,1)  -  R  .  C0S(.5  •  ANGL) 

2  S0(I,1)  -  R  •  SIN( . 5  »  ANGL) 

Cl  III! IF  OUTPUT  OF  UNWRAPPED  COORDINATES  IS  DESIRED 
C  WRITE  (6, 1000) 

C  WRITE  (6,2000)  ( I , A0( 1 , 1 ) , S0( 1 , 1 ) ,  I  -  1  .  II) 

RETURN 

1000  FORMAT (IX, ’UNWRAPPED  COORDINATES  IN  THE  TRANSFORMED  PUNE’) 
2000  FORMAT(I5  .  2F20.8) 

END 

SUBROUTINE  TABINT(XP, YP.XSING. YSING.N) 

DIMENSION  XP( 161 ) , YP( 161 ) , S0( 161 ) , A0( 161 ) 

C I  I  1 1 ! SMOOTH  THE  AIRFOIL  SURFACE  BY  FINDING  ADDITIONAL  POINTS 
U  -  XP(1)  -  XSING 

V  -  YP(1)  -  YSING 
U  -  1  . 

V  -  0. 

ANGL  -  8.  •  ATAN(  1  .  ) 

DO  1  I  -  1.N 
XII  -  XP(I)  -  XSING 
Y 1 1  -  YP(I)  -  YSING 

ANGL  -  ANGL  +  ATAN2( (U»Y1 1-V»X1 1 ) , (U»X1 1+V»Y1 1 ) ) 

R  -  SQRT(X1 1  »«2  +  Y11  •*  2) 

U  -  XII 

V  -  Y11 
R  -  SQRT(R) 

A0 ( I )  ■  R  •  COS (ANGL  •  .5) 

1  S0(I)  -  R  •  SIN(ANGL  •  .5) 

OX  -(A0(N)-A0(1))/96. 

A0ST  -  A0(1) 

DO  3  I  ■  1  ,97 
XX  -  FLOAT ( 1-1 )  .  OX  +  A0ST 
CALL  TAINT(A0 , S0, XX , YY ,N, 3 . NCR , MON) 
xpm  -  XX  •  XX  -  YY  •  YY  +  XSING 
3  YP(I)  -  2.  •  XX  •  YY  +  YSING 
RETURN 
END 

SUBROUT  I NE  TA I NT ( XTAB . FT AS , X . FX , N , K , NER . MON) 

DIMENSION  XTAB(1),FTA8(1),T(10),C(10) 

NASA  -  AMES  SUBROUTINE  FOR  POLYNOMIAL  INTERPOUTION 
OF  A  TABUUTED  FUNCTION 

IF(M-K)  1,1,2 

1  NER  -  2 
RETURN 

2  IF(K-9)  3,3,1 

3  I F (MON)  4,4,5 
5  IF(MON-2)  6,7,4 

4  J  -  0 
NM1  >  N  -  1 
DO  8  I  -  1  ,  NM1 
IF(XTA8( I )  -  XTAB( 1+1 ) )  9,11,10 

11  NER  -  3 
RETURN 
9  J  *  J-1 
GO  TO  8 
10  J  -  J+1 
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8  CONTINUE 
MON  -  1 

IF(J)  12  .  6  ,  6 

12  MON  -  2 

7  00  13  I  -  1  ,  N 

IF(X  -  XTAB(I))  14.14.13 

14  J  -  I 

GO  TO  18 

13  CONTINUE 
GO  TO  15 

6  DO  16  I  -  1  .  N 

IF(X-XTA9(I))  16.17.17 

17  J  -  t 

GO  TO  18 
16  CONTINUE 

15  J  -  N 

18  J  -  J  -  (K+1 )  /  2 
IF(J)  19,19,20 

19  J  -  1 

20  M  •«  J  +  K 

IF(M  -  N)  21 ,21 .22 

22  J  -  J  -  1 
GO  TO  20 

21  KP1  -  K  +  1 
JSAVE  -  J 

26  DO  23  L  -  1 ,  KP1 
C( L)  -  X  -  XTAB(J) 

T(L)  -  FTAB(J) 

23  J  -  J+1 

DO  24  J  -  1.K 
I  ■  J+1 

25  T(I)  -  (C(J)*T(I)<(I)«T(J))/(C(J)-C(I)) 

I F( I-KP1 )  25.25.24 

24  CONTINUE 
FX  -  T(KP1 ) 

NER  -  1 
RETURN 

END 

SUBROUTINE  SING<N2.N,X,Z.XIE,YIE.TEA.TES,XSING.YSING.CHD) 


THIS  SUBROUTINE  COMPUTES  SINGULAR  POINT  LOCATIONS. 
DIMENSION  X(2)  .  Z(2) 


NU 

-  N2 

hi 

*  Hi  +  | 

NJ 

■  HI  -  1 

Tl 

-  KfNli 

Z1 

*  l  Ml) 

1 2 

-  *(U21 

72 

-  Z(NJ) 

-  <£M 

n 

-  I(H Z) 

m 

-  M  «■  »  -  XI 

•  » 

2 

DJI 

*  u  <*  i  -  zi 

t* 

2 

01 

-  «  -  XI 

D4 

-  r?  -  n 

05 

*  Hi  **  2  -  XI 

•  • 

2 

D0 

-  73  ■■  r  -  zi 

•  • 

2 
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D7  -  X3  -  XI 
08  »  Z3  -  Z1 

8  -  f07  •  (  01  +  02)  -  03«(05+06))/(2.  *(D7«D4-03»D8)) 

IF(ABS(D3).LT.ABS(D7))  GO  TO  18 
A  -  (01  +  02  -  2.  •  B  •  04)  /  (2.  *  03) 

GO  TO  20 

10  A  -  (D5  +  06  -  2.  *  B  «  08)  /  (  2.  •  D7) 

20  CONTINUE 

R  -  SORT((X2-A)*»  2  +  (Z2-B).*2) 

XLE  -  X(NU) 

YLE  -  Z(NU) 

CHO  -  X(1)  -  X(NU) 

A2  -  (Z  2  -Z(1))  /(X(2)  -  X( 1 ) ) 

A1  -  (Z(N)-Z(N-I ) )/(X(N)-X(N-1 )) 

TES  -  .5  •  (At  +  A2) 

TEA  -  A2  -  A1 

TEA  -  TEA  .  57.29578 

XSING  -  (A+XLE)  /2. 

YSING  -  (8+YLE)  /  2. 

RETURN 

END 


SUBROUTINE  AIRFOL( 1 1 , JJ . IT. IE) 

COM40N/GRID1/X(161 ,41 ) ,Z( 161 ,41 ) 

COMMON/YSYM/ISYM 

DIMENSION  S0( 1 61 ,4) , A0( 161 ,4)  ,Y0(41 ,4) ,XP( 161 ) , YP( 161 
1E(161),F(161),B0(49) 


). 


DATA  (B0(  I ) , 1-1 ,32)/1 . ,  1.0414.1.0836,1 .1270,1.1715,1.2175. 1.2651 . 
11 .3145,1 .3659,1 .4199,1.4755.1 .5349.1.5973,1 .6636,1 .7342.1 .8099, 
21.8914,1 .9799,2.0764,2.1829.2.3012,2.4341 ,2.5653.2.7597,2.9646, 
33.2106.3.5141 ,3.9019.4.4219.5.1687.6.3632.8.6809/ 


III!  (COMPUTE  THE  COMPUTATIONAL  GRID  POINTS  BASED  ON  INPUT  AIRFOIL  SHAPE 
00  8  I  ■  1  ,32 
8  A0( 1,1)  .  B0(  I ) 

READ  (5.1) 

1  FORMAT(IX) 

READ  (5.2)  FNU.FNL.ZSYM 

2  FORMAT(JF10.0) 

ISYM  -  0 

IF(ZSYM.NE.0.)  ISYM  -  1 

11-157 

JJ  -  41 

IT  -  31 

IE  -  127 

IIP1  -  II  +  1 

IIM1  -  II  -  1 
1 1 JJ  -  1 1  *  JJ 
IIJJ2  -  II  .  (  JJ-2) 

ILE  -  (IT  +  IE  )  /  2 

ISTP  -  0 

NN  -  5 

NRF  -  0 

NOTAPE  -  1 


PI 

-  4. 

•  ATAN( 

NU 

-  FNU 

NL 

-  FNL 

N 

-  NU  + 

NL  -  1 

READ(5. 1 ) 

READ  (5,333)  (XP(I).YP(I) , I  -  NL  ,  N) 
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333  FORMAT(2F10.0) 

9994  FORMAT (F20. 8) 

L  -  N  +  1 

!F(ZSYM  .GT.  0.)  GO  TO  9995 

L  -  NL  +  1 

REAO(5,1) 

READ  (5,333)  (XP(L-t ) . YP(L-I ) , 1-1 ,NL) 
GO  TO  9996 

9995  K1  -  L 

DO  16  I  -  NL  .  N 
K  -  K1  -  I 


XP(K 

YP(K 


)  “  XP(I) 

)  -  '  VP(I) 


16  CONTINUE 
9996  SCALE  -  1.  /(XP( 1 )-XP(NL) ) 
XX  -  XP(  ' 


S(NL) 

9(NL) 


YY  -  YP( 

DO  9997  I  -  l‘  ,  N 
XP(I)  -  XPfl)  •  SCALE 
9997  YP( I )  ■  YP( I )  .  SCALE 

CALL  SING(NU,N, XP, YP , XLE , ZLE, TEA, TES , XSING, YSING.CHO) 

CALL  TABINT(XP,YP,XS1NG,YSINI!.N) 

NBOOY  -  IE  +  1  -  IT 
DO  6791  I  -  1  ,  NBOOY 
L  -  I  -  1 
EfIT+L)  -  XP(I) 

6791  F(IT+L)  -  YP( I ) 

IEP1  -  IE  +  1 
SLOPT  -  TES  •  .75 
DO  438  I  -  IEP1  .  II 
II  -  I  +1  -  IE 
E( I )  -  A0( 11,1) 

OXI  -  t.  /  48. 

0  -  4.  /  3.  •  (E( I )  -  .25) 

F(I)  -  F( IE)  +  SLOPT  .  ALOG(D)  /  D 
L  -  IIP1  -  I 
E(L)  -  Efl) 

F(L)  -  F( IT)  +  SLOPT  .  ALOG(D)/D 
WRITE  (6.439) 

FORMAT(2X.3H  I . 19X, 1HX, 19X. 1HY) 

WRITE  (6.37)  (I.E(I).F(I).I  -  1  .  II) 

CALL  WRAP ( 1 1 . J  J , XS I NG . YS I NG .  E , F , S0 , A0 , Y0) 

IMAP  GRID  BACK  TO  PHYSICAL  PLANE  AND  SHIFT  TO  QUARTER  CHORD 
DO  10  J  -  2  .  JJ 


438 


439 


Cl  It  I 


C 

c 

c 


DO  10  I  -  1  .  II 

1X(^J-’)  •  A0(  1 , 1  )««2  -  (S0(  1 , 1  )+Y0( J  ,  1 )  )»»2 

10  Z(I.j-l)  -  2.  *  A0( 1,1)  .  (S0( I . 1 )+Y0( J , 1 ) ) 
JJ  ■  JJ  —  1 
RETURN 

37  FORMAT( I5.2F20.8) 

END 

SUBROUTINE  CLUSTR(DS) 

C0M40N/GR ID1/X(161 ,41) ,Z(181 ,41) 
COMJON/DCRID/DT. IMAX.KMAX, ITEL, ITEU 
DIMENSION  S(41) ,XP(41),YP(41) ,R(41 ) 


THIS  SUBROUTINE  CLUSTERS  A  GIVEN  X 
THE  USER-SPECIFIED  DISTANCE  DNMIN 


Z  GRID  SUCH  THAT  THE  FIRST  POINT 


CUM  (COMPUTE  THE  OLD  STRETCHING 


S  AT 
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DO  100  I  -  1  ,  I MAX 
S(1)  -  0. 

XP(1)  -  X(I.I) 

YP(1)  -  2(1.1) 

DO  10  K  -  2  .  KMAX 
XP(K)  -  X( I , K) 

YP(K)  -  Z(l,K) 

10  S(K)  -  SQRT(  (XP(K)-XP(K-I  ))*»2+(YP(K)— YP(K-1 ) )»»2) 
1+S(K-1 ) 

SUMOX  -  S(KMAX) 

CALL  STRTCH(SUMDX,DS,F1 , KMAX , FACTOR) 

C  WRITE  (6,200)  I, FACTOR 

R(1)  -  0. 

OR  -  DS 

DO  20  K  -  2  .  KMAX 
R(K)  -  R(K-1)  +  DR 
DR  -  DR  •  FACTOR 
20  CONTINUE 

RUST  -  1  .  /  R(KMAX) 

DO  30  K  -  2  ,  KMAX 

R1  -  R(K)  .  RLAST  .  S(KMAX) 

Cl  I  I  I  (REDISTRIBUTE  THE  CONSTANT-ETA  LINES 

CALL  TAINT(S,XP,R1  ,XP1 .KMAX . 3 , NER ,MON) 

X(I.K)  -  XP1 

CALL  TAINT(S,YP,R1 . YP1 , KMAX , 3 , NEP , MON) 

Z(I.K)  -  YP1 
30  CONTINUE 
100  CONTINUE 
C  WRITE  (6,115) 

DO  110  I  -  1  ,  IMAX 
OX  -  X( ! , 2)  -  X( I , 1 ) 

DY  -  Z(  ,2)  -  Z ( I . 1 ) 

DN  -  SCRT(DX.OX+DY»OY) 

C  WRITc(6, 120)  I  .  OX  ,  DY  ,  DN 

110  CONTINUE 
RETURN 

115  FORMAT (5X.6HNORMAL,  1X.8HDISTANCE.3H  AT,4H  THE.5H  WALL,/ 
1.5H  I.ax^HOX.SX^HDY.SX^HDN,//) 

120  FORMAT (I5.3F10.5) 

200  FORMAT ( 15, F10 . 5) 

END 

SUBROUTINE  STRTCH( SUMOX, DX1  ,F1  ,N1  ,R) 


THIS  SUBROUTINE  DETERMINES  A  GEOMETRIC 
PROGRESSION  OF  GRID  SPACING  BETWEEN  1  '  ND  N1  SUH  THAT 
SUM8DX)  EQUALS  SLiCX.  THE  RATIO  BETWE'N  SUCCESSIVE 
SPACINGS  IS  R. 

N  -  HI  -  1 
R  -  1.5 


El  -  1.E-04 
E2  -  1 . E-04 
DO  10  L  -  1 .  50 

F-  (R— 1 )  •  SUMOX  -  DX1«(R**N-1) 

FP  -  SUMOX  -  DX1  .  FLOAT(N)  .  R  ••  (N-1) 
RITER  -  R  -  F /  FP  ’ 

IFfl  E-02.LT. RITER. AND. RITER.LT. 1.)  RITER  -  1 
im..LT  RITER.ANO.RITER.LT.  100. )  RITER-. 01 
IF(AB?(f  RITER). LT.  R«E1)  GO  TO  1 
n  -  R>  i  .R 
10  CONTINUE 
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R  -  1.0601 

0X1  -  DZTOT/FLOAT(N1-1) 

RETURN 
1  R-  RITER 
RETURN 
ENO 

SUBROUTINE  EDDY(DALFA) 

COMMON/FLOW/Q1  (161,41), 02(161, 41), 03(161  .41 ) ,Q4( 101 , 41 ) 
C0M40N/MUTUR/CMU (161,41) 

C0M40N/SK I NCF/CF( 161) 

COMMON/DCR I D/DT , I MAX , KMAX , I T EL .  I T  EU 

COMMON/P AR/GAVWA .  REYREF ,  ALFA .  ALFA1 .  REDFRE ,  AMINF .  ALFA  I 

COMMON/GR I 0 1 /X ( 101 , 41 ) ,Z( 181 ,41) 

DIMENSION  T IN(41 ) , TOUT (41 ) , Y(41 ) 

INITIALIZE  VISCOSITY  EVERYWHERE 
FACT1  -  OT  •  AMINF  /  REYREF 
CMUMAX  -  100.  •  FACT1  /  DT 
DO  1  K  -  1  ,  KMAX 
DO  1  I  -  1  ,  IMAX 
1  CMU(I.K)  -  FACT1 

THIS  SUBROUTINE  COMPUTES  THE  EDDY  VISCOSITY  BASED  ON  THE 
BA LOW IN- LOMAX  TWO  LAYER  MOOEL 


C 


DO  100  I  -  2  ,  IMAX  -  1 
UDIF  -  0. 

FMAX  -0.1 E-06 
YMAX  -  . 1 E-06 
FYMAX  -  0. 


Y(1)  -  0. 

UWALL  -  0. 

IF( I .GT . ITEU.OR. I . LE  ITEL)UWALL  -  SQRT(Q2(I . 1 )»*2+Q3( I , 1 )»*2)/ 

101(1,1) 

COMPUTE  TAU  AT  THE  WALL 


UET 

m 

1 . *(02( I ,2)/Q1 (1,2 
1 . •(03( I ,2)/Q1 (1,2 

VET 

m 

XXI 

m 

X(I+1.1)  -  X( 1-1 . 1 

ZXI 

m 

Z(I+1.1)  -  Z< 1-1 . 1 

XET 

m 

4.  .  X( I , 2)  -  3.  • 

ZET 

m 

4.  .  Z( I ,2)  -  3.  • 

XXI 

m 

.5  •  XXI 

ZXI 

m 

.5  •  ZXI 

XET 

m 

.5  •  XET 

ZET 

m 

.5  •  ZET 

-  Q2(I.1)/Q1(I.1)) 

-  Q3( I , 1 )/Q1 (1.1); 


YAC  -  1.  /  (XXI  •  ZET  -  ZXI  .  XET) 

OMEGA  -  (UET  .  XXI  -  VET  •  ZXI  )  •  YAC 
TWALL  -  AMINF  •  OMEGA  /  REYREF 
CF(I)  -  2.  •  TWALL  /  (AMINF««2) 

FACT  -  SQRT(Q1 (1,1)  •  ABS(TWALL))*REYREF/(20. *AMINF) 
DO  10  K  -  2  ,  KMAX— 1 


UXI  -  (Q2( 1+1 ,K)/Q1 ( 1+1 ,K)  -  02(1-1 ,K)/Q1 ( 1-1 ,K)) 
VXI  -  (Q3( 1+1  ,K)/Q1 (1  +  1 ,K)-Q3(  1-1 ,K)/Q1 ( 1-1 ,K)) 

UET  -  (Q2( I ,K+1 )/Q1 ( I ,K+1 )— 02 ( I  ,K-1 )/Q1 ( I ,K-1 )) 

VET  -  (Q3( I ,K+1 )/01 ( I ,K+1 )-Q3( I ,K-1 )/Q1 ( I ,K— 1 )) 

XXI  -  X(I+1 ,K)  -  X( I— 1 ,K) 

ZXI  -  Z(I+1 ,K)  -  Z(I-1 ,K 
XET  -  x( I . K+1 )  -  X(I.K-l) 

ZET  -  Z(I ,K+1)  -  Z(I ,K-1) 

YAC  -  1.  /  (XXI  •  ZET  -  ZXI  •  XET) 

OMEGA  -  ABS ( UET • XX I +VET • ZX I -UX I  »X ET-VX I • ZET )  ♦  YAC 
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UDIF  -  SORT (Q2( I ,K)**2+Q3( I ,K)*»2)/Q1 ( I ,K)  -  UWALL 
IF(ABS(UOIF) .GT.UOIFMAX)  UOIFMAX  -  ABS(UOIF) 

Y(K)  -  SORT ( (X( I ,K)-X( I ,K-1 ) )««2+(Z( I ,K)-Z( I ,K-1 ) )*»2)+Y(K-1 ) 

F  -  YfK)  .  OMEGA 
IF((Y(K)»FACT) .GT.20. )  GO  TO  31 

IF(I.GT.ITEL.ANO.I.LT.ITEU)  F  -  F  *  (1.  -  EXP(-Y(K)»FACT) ) 

31  CONTINUE 

MOOIFIED  TURBULENCE  MOOEL  APPLY  FOR  SPECIFIED  RANGE  OF  ANGLES  WHERE 
FY  IS  USED  T0  FI NO  THE  SECONO  PEAK  VALUE  OF  F  FUNCTION 

IF(ALFA. LT.ALFAI . ANO.DALFA.GE.0.  )  THEN 
FY  -  F  .  Y(K) 

IF(FY.GT.FYMAX)  THEN 
FYVtAX  -  FY 
FMAX  -  F 
rMAX  -  Y(K) 

END  IF 
END  IF 

IF(ALFA.GE.ALFAI . OR.DALFA. LT.0. )  THEN 
IF(F.GT.FMAX)  THEN 
FMAX  -  F 
YMAX  -  Y(K) 

END  IF 
END  IF 

FCT  -  Y(K)  •  FACT 
IF(FCT. GT.20.)  FCT  -  20. 

FCT  -  ABS(FCT) 

EL  -  .4  •  Y(K)  •  (1.  -  EXP(-FCT) ) 

TIN(K)  -  QI(I.K)  •  EL  •  EL  *  OMEGA 
TIN(K)  -  ABS(TIN(K) ) 

10  CONTINUE 
KSWTCH  -  0 
FWAKE  -  YMAX  •  FMAX 
FI  -  0.25  •  YMAX  .  UDIF  ..2  /  FMAX 
I F(F1 . LT . FWAKE)  FWAKE  -  FI 
DO  20  K  -  2  .  KMAX  -  1 
FKLEB  -  0. 

I F( ABS( Y(K)/YMAX) . LT . 1 . E+04)  THEN 
FKLEB  -  1.  /  (1.  +  5.5  •  (0.3  ♦  Y(K)/YMAX)  ••  6) 

END  IF 

TOUTfK)  -  .0168  •  1.6  «  QI(I.K)  •  FWAKE  *  FKLEB 
TOUT(K)  -  ABS(TOUT(K) ) 

IFfKSWTCH.NE.0)  GO  TO  20 
IF(TIN(K).GT. TOUT(K) )  KSWTCH  -  K  -  1 
20  CONTINUE 

CM  !  !  (TOTAL  VISCOSITY  IS  SUM  OF  LAMINAR  AND  INNER/OUTER  LAYER  AS  APPROPRIATE 
DO  30  K  -  2  ,  KMAX  -  1 
IF(K. LE. KSWTCH)  THEN 

OAJ(I.K)  -  DT  •  (AMINF/REYREF  +  ABS(TIN(K))) 

ELSE 

CMU(I.K)  -  DT  .  (AMINF  /  REYREF  +  ABS(TOUT(K) ) ) 

END  IF 
30  CONTINUE 
100  CONTINUE 
RETURN 
ENO 

SUBROUTINE  RES  I (OMEGA) 

C0U4ON/PERTR/DQ1(161 .41 ) ,DQ2( 161 , 41 ) ,D03( 161 ,41 ) ,004(161 .41) 
C0M40N/GRID1/X(161 . 41 ) , Z( 1 61 ,41) 
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C0MM0N/DGR1D/DT , IMAX, KMAX.  ITEL.ITEU 

C0M40N/FLOW/Q1 (161 ,41 ) ,02(161 ,41 )  ,G3(  161 . 41 ) ,Q4( 161 .41 ) 
C0M40N/PAR/GAMAA ,  REYREF .  ALFA .  ALFA  1 ,  REOFRE .  AMI  NF ,  ALFA  I 
DIMENSION  RHS( 161 , 4) 

XTAU(I.K)  -  OMEGA  •  Z(I.K) 

YTAU(I.K)  -  -  OMEGA  •  X(I,K) 

THIS  SUBROUTINE  COMPUTES  THE  RESIDUAL  ON  THE  RIGHT  HANO 
SIDE  ARISING  FROM  THE  EULER-  PART  OF  THE  GOVERNING  EQUATIONS 

FLUX  TERMS  WITHIN  THE  XI-  DERIVATIVE 
DO  100  K  -  2  ,  KMAX  -  1 
DO  10  I  -  1  .  IMAX 

UCON  -  (Q2( I ,K)/Q1 ( I ,K)  )  «  (Z( I ,K+1 )-Z( I , K— 1 ) ) 

1  -  (Q3( I ,K)/Q1 ( I ,K) )  .  (X{ I ,K+1 )— X ( I , K-1 ) ) 

UCON  -  0.25  •  DT  •  UCON 

XIT  -  -  XTAU(I.K)  • ( Z ( I , K+ 1 ) — Z ( I , K— 1 ) ) 

1  +  YTAU( I , K)  •  (X(I,K+1)  -  X(l.K-l)) 

XIT  -  XIT  •  OT  •  0.25 
UCON  -  UCON  +  XIT 
RHS(I.I)  -  UCON  •  QI(I.K) 

R  -  1 .  /  QI(I.K) 

P  -  (GAMUA-1.)  .  (Q4( I ,K)  -  .5  •  R  • (02(1 ,K)»»2+ 

1  Q3( I ,K)*»2)) 

RHS( I ,2)  -  02(1, K)  .  UCON  +  P  *  DT  •  0.25  *  (Z(I,K+1)  -  Z( I .K-1 ) ) 
RHSU.3)  -  03(1.  K)  •  UCON  -  P  •  DT  •  0.25  •  (X( I ,K+1 )-X( I ,K-1 ) ) 
RHS( I . 4)  -  UCON  .  (Q4( I ,K)+P)  -  XIT  •  P 

10  CONTINUE 

DO  11  I  -  2  ,  IMAX  -  1 

DQI(I.K)  -  DQI(I.K)  -  RHS( 1+1 , 1 )  +  RHS(I-1,1 
DQ2U.K)  -  DQ2  ( I .  K )  -  RHS(I+1,2)  +  RHS(I-1.2 
DQ3(I.K)  -  003 ( I .  K  )  -  RHS(I+1,3)  +  RHS(I-1.3 

11  D04( I ,K)  -  DQ4( I  ,K)  -  RHS( 1+1 .4)  +  RHS(I-1,4 
100  CONTINUE 

FLUX  TERMS  WITHIN  THE  ETA-  DERIVATIVE 

DO  200  I  -  2  ,  IMAX  -  1 
DO  20  K  -  1  ,  KMAX 
VCON  -  (Q2(I.K)/Q1(I,K) 

1  +(03 ( I ,K)/Q1 ( I ,K) 

VCON  -  VCON  •  0.25  •  DT 
ZTAT  -  -XTAU(I.K)  .  (Z(  1-1  ,K)-Z( 1+1 ,K) )  -  YTAU(I,K). 

1  (X( 1+1 ,K)-X(I-1 ,K)) 

ETAT  -  ETAT  •  0.25  •  DT 
VCON  -  VCON  +  ETAT 
RHS(K.I)  -  VCON  «  01 ( I ,K) 

P  -  (GAM4A— 1 . )  .  (Q4f I . K I  -  0.5  *(Q2( I ,K)»»2+Q3( I ,K)»»2)/Q1 (I , K ) ) 
RHSfK.2)  -  VCON  .  Q2( I ,K )  +P  •  DT  .  .  25  •  (Z( 1-1 ,K)-Z( 1+1 ,K) ) 
RHS(K,3 1  -  VCON  •  03(1, K)  +  P  •  DT  .  .  25  •  (X(I+1,K)  -  X(I-I.K)) 
RHS(K ,4)  -  VCON  •  (04(1, K)+P)  -  ETAT  •  P 

20  CONTINUE 

DO  21  K  -  2  .  KMAX  -  1 

DQI(I.K)  -  DQI(I.K)  -  RHS(K+1 , 1 )  +  RHS(K-1 , 1 
002(1. K)  -  002(1. K)  -  RHSiK+1 ,2)  +  RHSCK-1,2 
003(1, K)  -  DQ3( I ,K )  -  RHS(K+1,3)  +  RHSCK-1,3 

21  D04(l.K)  -  D04(l,K)  -  RHS(K+1,4)  +  RHS(K-1,4 
200  CONTINUE 

RETURN 
END 

SUBROUT  I NE  ROTGR 1 0 ( X , Z , I MAX , KMAX , DALFA ) 
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C  ROTATE  GRID  IN  THE  CLOCKWISE  DIRECTION  BY  AN  AMOUNT  DALFA 
DIMENSION  X(161 ,41 ) ,Z(161 ,41 ) 

CA  -  COS (DALFA) 

SA  -  -  SIN(DALFA) 

00  20  K  -  1  .  KMAX 

DO  20  I  -  1  .  IMAX 

XI  -  X ( I , K ) 

21  -  2(1 , K) 

X(I ,K)  -  XI  •  CA  -  Z1  •  SA 

20  2(1 .K)  -  Z1  •  CA  +  XI  •  SA 

RETURN 
END 

SUBROUTINE  CPPLOT (11,12, FMACH , X , Y , CP ) 

C 

C  THIS  SUBROUTINE  PLOTS  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE 
C 

DIMENSION  K0DE(4) ,LINE(90),X(161),Y(161),CP(161) 

DATA  K0DE/1H  . 1H+, 1HI , 1H*/ 

WRITE  (6.2) 

2  FORMAT(50H0PLOT  OF  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE/ 

1  10H0  X  .  10H  X/C  .  1 0H  CPL  ,  1 0H  CPU  ) 

CP0  -  (1.  +  .2  •  FMACH  *«2)  ••  3.5  -  1. 

CP0  -  CP0  /  (  .7  «  FMACH  »»2) 

K0  -  30.  •  CP0  +  4.5 
1MIN  -  ( 1 2-1 1 )/2  +  II 
I  LOW  -  2  *  IMIN 
CHD-X(II)  -  X(IMIN) 

DO  12  I  -  1  ,90 
12  LINE(I)  -  KOOE( 1 ) 

DO  34  I  -  IMIN  .  12 
K  -  30.  *  (CP0  -  CP( I ) )  +  4.5 
K1  -  30.  •  (CP0  -CP ( I  LOW- 1 ) )  +  4.5 
IF(K.LT.I)  K  -  1 
IF(KI.LT.I)  K1  -  1 
I F(K .GT . 90)  K  -  90 
I F(K1  GT.  90)  K1  -  90 
LINE(K0)  -  K00E(3) 

LINE(K)  -  K00E(2) 

LINE(KI)  -  K00E( 4) 

XOC  -  (X(I)  -  X( IMIN) )  /  CHO 

WRITE  (6,610)  Xfl), XOC. CP(ILOW-I),CP(I),  LINE 

LINE(KI)  -  KOOE( 1 ) 

34  LINE(K)  -  KOOE( 1 ) 

RETURN 

610  FORMAT(4F10.4,90A1 ) 

ENO 

/EOF 

NACA  0012  AIRFOIL.  RUN: 

'5®  -0025  5.  15.000  10.00  0.00  0.151 

NO.  OF  STEPS 
4500. 

REYNOLDS  NUMBER  IN  MILLIONS,  OISTANCE  OF  FIRST  POINT  OFF  THE  WALL 
345  . 00005 

TSTART 

FULOUT 
-1  .0 

RESTART, PITCH, PLUNGE. OUTPUT  OPTIONS  SPECIFIED  IN  THE  NEXT  CARD 
TRUE  TRUE  FALSE 

TNU  FNL  FSYM 


.5000 


111 


33. 

33. 

X 

Y 

e. 

0. 

0.0050 

.01153 

.0125 

.01894 

.0250 

.02615 

.0500 

.03555 

.0750 

.04200 

.1000 

. 04683 

.1500 

. 05345 

.2000 

.05737 

.2500 

.05941 

.2600 

.05962 

.2700 

.05978 

2800 

.05992 

2900 

.05999 

.3000 

.06000 

3100 

05999 

.3200 

.05992 

.3300 

. 05980 

.3400 

. 05965 

.3500 

. 05947 

.4000 

.05803 

.4500 

.05581 

.5000 

.05294 

.5500 

. 04952 

.6000 

. 04563 

.6500 

.04137 

.7000 

.03664 

.7500 

.03161 

.8000 

.02623 

.8500 

. 02055 

.9000 

.01448 

.9500 

.00807 

1 . 0000 

.00126 
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APPENDIX  B 

NOTES  ON  USE  OF  THE  NAVIER-STOKES  SOLVER 
1.  JOB  CARDS 

The  JCL  options  are  selected  by  removing  the  "comment"  designator  ("*.")  at  the 
beginning  of  the  applicable  lines.  The  options  available  are: 

•  Save  the  current  solution.  Values  of  TIME,  Ql,  Q2,  Q3,  and  Q4  are  saved 
(logical  unit  08)  for  subsequent  restart.  Activate  the  two  lines  referencing 
"NEWS  LX". 

•  Create  pressure  coefficient  data  file.  Data  is  accumulated  through  the  runs  to 
be  accessed  and  read  by  a  separate  program.  Activate  the  "OLDCP"  statements 
to  access  and  add  to  previously  stored  data.  Activate  "XEWCP"  statements  to 
store  current  cumulative  data.  A  new  version  is  created  each  time,  so  the  files 
must  be  purged  periodically. 

•  Start  from  a  stored  solution.  The  above  values  are  read  (logical  unit  07)  and 
iteration  continued  from  that  point.  Activate  the  two  lines  referencing 
"OLDSLX". 

•  Creat  PLOT3D  files.  Conversion  from  Cray  to  VAX  binary  is  handled  and 
properly  formatted  "Q"  and  "XYZ"  files  are  created  for  use  with  the  PLOT3D 
graphics  program.  Activate  the  Q  and  XYZ  "ASSIGN"  statements  and  the 
"SEXDVAX"  and  "ACCESS"  statements  before  the  "EXIT"  card. 

•  Create  "troubleshooting"  PLOT3D  files.  If  the  solution  "blows  up"  and  the 
appropriate  "WRITE"  statements  are  not  commented  out  in  the  main  program, 
Q  and  XYZ  files  are  created  to  investigate  the  status  of  the  solution  at  the  last 
"successful"  iteration.  Activate  the  "ASSIGN"  statements  and  the  "ACCESS" 
and  "SEXDVAX"  statements  after  the  "EXIT"  card. 

•  Use  job  chaining.  The  FETCH,  REWIND,  and  SUBMIT  statements  are  used 
to  call  up  another  program  when  current  run  is  completed. 

In  addition,  if  it  is  desired  to  save  more  than  one  solution,  the  PDX  (5  digits) 
and  ID  may  be  changed.  Old  data  sets  must  be  purged  from  Cray.  This  is 
accomplished  by  placing  an  AUDIT.'  card  after  the  account  card  when  running  a 
program  or  by  running  "AUDIT.JCL"  to  obtain  a  listing  of  current  data  sets.  Then, 
on  VAX,  run  SKILLJCL  to  create  a  JCL  to  delete  Cray  PDN's.  It  is  then  only 
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necessary  to  precede  this  JCL  by  a  JOB  card  and  ACCOUNT  card  and  run  the 
program  as  usual.  Consult  the  ACF  Cray  Users'  Guide  for  detailed  information  on  job 
control  cards. 

2.  MAIN  PROGRAM 

Certain  changes  may  be  made  within  the  main  program.  These  changes  affect 
the  program  execution  or  change  the  output. 

•  The  frequency  of  steady-state  output  may  be  changed  by  varying  the  value  in 
the  first  statement  after  the  comment  "FOR  STEADY  STATE  OUTPUT  USE 
THE  FOLLOWING". 

•  The  interval  for  exiting  the  program  (in  order  to  save  a  solution  and  generate 
PLOT3D  files)  may  be  changed  under  "MULTIPLY  NCPOUT...".  For 
example:  NEXIT  =  24  *  NCPOUT  means  the  program  will  exit  automatically 
four  times  a  cycle,  or  every  24  times  the  normal  printed  output  is  generated. 

•  Velocity  profile  information  may  be  output  if  desired  (as  when  a  permanent 
record  of  a  converged  solution  is  desired)  by  activating  the  WRITE  statement 
jus.  before  the  CONTINUE  statement  numbered  4000.  This  outputs  the 
indices  of  the  X  and  Z  coordinates  of  each  grid  point  with  the  corresponding 
values  of  pu,  pv,  eddy  viscosity,  total  velocity  (Jn2  +  v2),  and  distance  normal 
to  the  wall. 

•  At  the  beginning  of  subroutine  SLPS,  the  value  of  the  implicit  damping 
coefficient  may  be  adjusted  by  changing  the  multiple  of  WW. 

•  At  the  end  of  subroutine  SLPS,  the  frequency  with  which  the  residuals  are 
output  may  be  set  by  changing  the  value  under  "SELECT  THE  INTERVAL...". 
At  least  every  ten  steps  is  recommended. 

•  Other  output  values  may  be  turned  on  and  off  as  well.  "UNWRAPPED 
COORDINATES  IN  THE  TRANSFORMED  PLANE"  is  in  subroutine 
"WRAP".  Airfoil  coordinates  are  in  "AIRFOL".  Grid  spacing  and  "NORMAL 
DISTANCE  AT  THE  WALL"  are  in  "CLUSTR". 

3.  DATA  CARDS 

Most  of  the  "tuning"  of  the  program  is  done  with  the  data  cards: 

•  1st  Line  -  Name  of  airfoil  and  other  identifying  information.  Eighty  characters 
available. 
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•  2nd  Line  -  (1)  and  (2)  The  first  two  values,  [MAX  and  KMAX  (format  15),  are 
the  number  of  X  and  Z  coordinate  locations  to  be  used.  These  remain  at  159 
and  41  for  the  161  *41  C-grid  used  at  present.  The  remainder  of  the  line 
contains  seven  values  in  format  F10.0.  (3)  DT:  size  of  the  time  step.  This  is 
automatically  set  to  1.0  within  the  program  when  the  reduced  frequency  is  less 
than  or  equal  to  0.001.  In  this  case  the  program  uses  the  local  time-step  option 
(relaxation).  For  dynamic  stall  DT  =  0.005  is  recommended,  at  angles  below 
"0  degrees,  smaller  at  higher  angles.  (4)  WW:  explicit  artificial  viscosity  term. 
Normally  2-5,  with  about  2-3  recommended  for  static  cases.  5  for  dynamic  stall. 
Higher  numbers  have  greater  effect  on  solution.  (5)  ALFA:  mean  angle  of 
attack.  (6)  ALFA1:  amplitude  of  oscillation.  (7)  ALFAI:  angle  below  which 
a  modified  turbulence  model  is  used  (upstroke  only)  to  compute  eddy  viscosity 
(Baldwin- Lomax  model).  Normally  set  to  19  degrees  for  dynamic  stall.  May 
affect  stability  of  solution.  (8)  REDFRE:  reduced  frequency.  (9)  AMINF: 
Mach  number. 

•  3rd  Line  -  "NO.  OF  STEPS"-not  read. 

•  4th  Line  -  FNSTP:  number  of  time  steps  to  be  done  on  the  present  run  (format 
F10.4). 

•  5th  Line  -  "REYNOLDS  NUMBER..."-not  read. 

•  6th  Line  -  (1)  REYREF:  the  Reynolds  number  in  millions  (format  F10.4).  A 
negative  value  means  inviscid  flow.  (2)  DNMIN:  distance  of  the  first  point  off 
the  wall  (format  F10.4).  For  Reynolds  numbers  up  to  3  million  0.00005  should 
be  used  normally.  Stability  of  the  solution  may  be  improved  by  increusing  this 
value  in  some  cases  (ie.,  high  AOA  steady  angle  of  attack). 

•  7th  Line  -  "TSTART"— not  read. 

•  8th  Line  -  TSTART  (format  F10.4):  time  calculations  have  been  advanced  up 

to  the  previous  run.  When  a  negative  value  is  used,  TSTART  is  read  from 

logical  unit  08  (see  JCL  comments).  Then  normally  use  0.0  for  initial  runs  and 
— 1.0  for  restarts.  Must  use  0.0  for  first  dynamic  run  from  converged  steady- 
state  solution. 

•  9th  Line  -  "FULOUT"— not  read. 

•  10th  Line  -  FULOUT  (format  F10.4):  —1.0  means  no  plotting  files  will  be 

generated.  Set  0.0  to  begin  full  output,  then  set  1.0  to  continue.  When  using 

full  output,  the  appropriate  job  cards  must  be  activated. 

•  1 1th  Line  -  "RESTART,  PITCH,. .."-not  read. 
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•  12th  Line  -  RSTRT,  PITCH.  PLUNGE  (format  L5):  If  RSTRT  is  set.  stored 
values  of  TIME,  Ql,  Q2.  Q3,  and  Q4  are  read  to  continue  iteration.  PITCH  is 
set  for  dynamic  stall  and  indicates  change  in  angle  of  attack.  PLUNGE  will  not 
be  set  for  present  purposes.  It  indicates  up  and  down  motion  of  the  airfoil. 

•  13th  Line  -  The  remaining  lines  define  airfoil  geometry  and  for  the  present  are 
set  for  the  NACA0012  airfoil. 

4.  ADDITIONAL  NOTES 

For  high  angle  of  attack  (separated  flow),  or  if  otherwise  desired  for  time 
accurate  steady-state  calculations,  use  reduced  frequency  of  0.002  and  ALFA1 
(amplitude)  =  0.  Set  PITCH  =  FALSE  and  DT  =  0.00\  as  for  dynamic  stall.  For 
stability  WW  should  be  set  to  5  (or  even  higher,  up  to  lO)  and  ALFAI  at  or  below  the 
minimum  angle  to  use  the  original  Baldwin- Lomax  turbulence  model.  The  distance  of 
the  first  point  off  the  wall  may  also  be  increased  (~  0.0001),  which  has  the  effect  of 
coarsening  the  grid  for  the  troublesome  fine-mesh  leading-edge  area. 

When  doing  dynamic  stall  simulations  where  the  AOA  goes  to  the  20-25  degrees 
range,  use  the  DT  =  0.005  at  lower  angles  for  computation  time,  but  reduce  this  value 
when  restarting  going  into  the  high  angle  portion  of  the  cycle. 

The  values  of  DRMAX,  DUMAX,  etc.  may  be  useful  when  a  solution  "blows 
up".  The  IR  and  KR  values  give  the  indices  of  the  X  -  Z  location  on  the  grid  where 
density  changed  the  fastest  (IR  =  80  is  the  leading  edge).  Normally,  this  will  be  near 
the  leading  edge.  DT  should  then  be  reduced. 

Reynolds  numbers  of  106  to  10  x  106  should  show  no  effect  on  sensitivity  of 
solution  for  distance  of  the  first  point  off  the  wall  of  .00005,  since  this  places  several 
points  in  the  boundary  layer. 

Mach  numbers  of  .1-.5  should  not  alter  calculation  time  significantly,  but 
perhaps  twice  as  many  iterations  may  be  required  for  convergence  in  the  .5-. 8  range. 
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